Key Results
Prep
# Load packages, source functions, set plot parameters, get behavioral and DIC dfs
packages_to_load <- c(
"dplyr",
"tidyverse",
"latex2exp",
"purrr",
"tidyr",
"rlang",
"patchwork",
"haven",
"RWiener",
"psych",
"data.table",
"lme4",
"lmerTest",
"foreach",
"beset",
"faux"
)
sapply(packages_to_load, require, character.only=TRUE)## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: tidyverse## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ stringr 1.4.0
## ✓ tidyr 1.2.0 ✓ forcats 0.5.1
## ✓ readr 2.1.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()## Loading required package: latex2exp## Loading required package: rlang##
## Attaching package: 'rlang'## The following objects are masked from 'package:purrr':
##
## %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
## flatten_lgl, flatten_raw, invoke, splice## Loading required package: patchwork## Loading required package: haven## Loading required package: RWiener## Loading required package: psych##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## Loading required package: data.table##
## Attaching package: 'data.table'## The following object is masked from 'package:rlang':
##
## :=## The following object is masked from 'package:purrr':
##
## transpose## The following objects are masked from 'package:dplyr':
##
## between, first, last## Loading required package: lme4## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loading required package: lmerTest##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step## Loading required package: foreach##
## Attaching package: 'foreach'## The following objects are masked from 'package:purrr':
##
## accumulate, when## Loading required package: beset##
## Attaching package: 'beset'## The following object is masked from 'package:psych':
##
## r2d## Loading required package: faux##
## ************
## Welcome to faux. For support and examples visit:
## https://debruine.github.io/faux/
## - Get and set global package options with: faux_options()
## ************##
## Attaching package: 'faux'## The following object is masked from 'package:rlang':
##
## %||%## The following object is masked from 'package:purrr':
##
## %||%## dplyr tidyverse latex2exp purrr tidyr rlang patchwork
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## haven RWiener psych data.table lme4 lmerTest foreach
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## beset faux
## TRUE TRUEsf <- function() sapply(fs, source)
fs <- c(
paste0('./Modules/', list.files('./Modules/'))
)
sf()## ./Modules/GenPPCs.R ./Modules/GenPPCs2.R ./Modules/Utils.R
## value ? ? ?
## visible FALSE FALSE FALSESetPlotPars() 8/2/22 - Replacing old paths with public access data throughout doc wherever data shareable (leaving out raw trace files because they’re too big)
# Get main SRET_CTL df..
#bx_df <- read.csv("./../../data/cleaned_files/s_bdf.csv")
bx_df <- read.csv("./Public_Data/s_bdf.csv")
# .. and df with some extra info to put in
#extra_df_full <- read.csv("./../../data/raw_files/tx_cond+ffmq+maas_from_KW_6.24.2020.csv")
extra_df_full <- read.csv("./Public_Data/tx_assignments.csv")4/24/22 - See non-public directory for demographics
Figure 1
DDM Vis
starting_point <- .05
non_dec_time <- c(.1)
threshold <- 0
decision_bounds <- c(max(non_dec_time), (max(non_dec_time)+1+threshold))
decision_time <- seq(decision_bounds[1], decision_bounds[2],
by=(decision_bounds[2]-decision_bounds[1])/100)
CreateDriftTrajectories <- function(upper_drift, lower_drift, n_trajects, label, sd_drift=.02) {
### Simulate n_trajects drift trajectories ###
dfs <- list()
for (nt in 1:n_trajects) {
this_upper_drift <- rnorm(1, upper_drift, .01)
this_lower_drift <- rnorm(1, lower_drift, .01)
traject <- c()
for (dec_ts in seq_along(decision_time)) {
y_incr <- rnorm(1, this_upper_drift, sd=sd_drift)
n_incr <- rnorm(1, this_lower_drift, sd=sd_drift)
incr_total <- sum(y_incr, n_incr)
if (dec_ts == 1) traject[dec_ts] <- starting_point
if (dec_ts > 1) traject[dec_ts] <- traject[dec_ts-1] + incr_total
}
traject[traject > 1] <- 1
traject[traject < -1] <- -1
dfs[[nt]] <- setNames(data.frame(traject, label, decision_time, nt),
c("trajectory", "drift_type", "step_index", "run"))
}
dfs %>% bind_rows()
}
drift_1 <- CreateDriftTrajectories(.06, -.01, 12, "high", sd_drift=.02)
drift_2 <- CreateDriftTrajectories(.06, -.03, 12, "medium", sd_drift=.02)
drift_3 <- CreateDriftTrajectories(.03, -.07, 12, "low", sd_drift=.02)
all_drifts <- rbind(drift_1, drift_2, drift_3)
all_drifts$drift_type <- factor(all_drifts$drift_type, levels=c("high", "medium", "low"))
rect <- data.frame(xmin=0, xmax=seq(non_dec_time, .2, .033), ymin=-1, ymax=1)
drift_sample <- ggplot(all_drifts, aes(step_index, trajectory, group=interaction(run, drift_type))) +
geom_line(size = 2, alpha=.5, aes(color=drift_type)) +
geom_hline(yintercept=c(-seq(1, 1.6, .1), seq(1, 1.6, .1)),
alpha=c(rev(seq(.4, 1, .1)), rev(seq(.4, 1, .1))), color="gray57", size=2) +
geom_segment(aes(x=0, xend=non_dec_time, y=.05, yend=.05),
inherit.aes = FALSE, color="black", size = 2) +
geom_segment(aes(x=0, xend=non_dec_time, y=.2, yend=.2),
color="gray57", size = 2) +
geom_segment(aes(x=0, xend=non_dec_time, y=-.10, yend=-.10),
inherit.aes = FALSE, color="gray57", size = 2) +
ga + ap + theme(legend.text = element_text(size = 18),
legend.title=element_text(size=20, face = "bold"),
legend.key.size = unit(1.5, 'lines')) +
scale_color_manual(name=" drift rate",
values = c("green", "orange", "red")) +
xlab("time") + ylab("") +
#ylim(0, -2) +
geom_rect(data = rect, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax),
alpha = .1, inherit.aes = FALSE) +
# xlim(0, xlim_ub) +
theme(axis.text=element_blank()) + theme(axis.ticks=element_blank()) #+
drift_sample
#ggsave("../paper/figs/drift_sample.png", drift_sample, width=11, height=6, dpi=700)Simulating densities — commented bc it takes a long time to run
# a <- 1.93
# t <- .9
# b <- .5
#
# #high_drift_samples <-
# diff_drift_samples <-
# rbind(
# data.frame(rwiener(5e4, alpha=a, tau=t, beta=b, delta =.7), type="high"),
# data.frame(rwiener(5e4, alpha=a, tau=t, beta=b, delta =.04), type="medium"),
# data.frame(rwiener(5e4, alpha=a, tau=t, beta=b, delta = -.7), type="low"))
#
# diff_drift_samples[diff_drift_samples$resp=="lower", "q"] <- -diff_drift_samples[diff_drift_samples$resp=="lower", "q"]
#
# #diff_drift_samples$resp=="lower"
# test_diff_plot <- ggplot(diff_drift_samples, aes(x=q, color=type)) +
# geom_density(size=3.5) + ga + ap + tol +
# scale_color_manual(
# values = c("green", "orange", "red")) +
# theme(axis.text = element_blank(), axis.ticks = element_blank()) + xlab("") + ylab("")
# test_diff_plot#ggsave("../paper/figs/diff_drift_dists_resize.png", test_diff_plot, width=35, height=3)Modeling
GetMapEsts <- function(s_df, format="long") {
### Returns dataframe of Maximum a Posteriori estimates for the variables in a dataframe.
# s_df must be a subject-level dataframe (ie. excluding group-level estimates) ###
# Handle entry of dfs other than subject-only traces
if (!all(grepl("subj.", names(s_df)))) stop("Must enter a trace data-frame with ONLY
subject-level vars")
# Posterior mean
col_means <- colMeans(s_df)
# Map est
if (format == "long") {
mdf <- foreach(i = 1:ncol(s_df)) %do% {
this_par <- s_df[, i]
mdf_row <- data.table(
"ID"=unlist(map(strsplit(names(s_df)[i], "subj."), 2)),
"var"=unlist(map(strsplit(names(s_df)[i], "subj."), 1)),
"map_est"=as.numeric(bayestestR::map_estimate(this_par))
)
mdf_row
} %>% bind_rows()
z_IDs <- grep("trans", mdf$ID)
if (length(z_IDs)) {
mdf[grepl("trans", mdf$ID), "ID"] <- unique(mdf$ID[!grepl("trans", mdf$ID)])
}
}
if (format == "short") {
# Put in short format with each variable as a column name
# First extract the unique vars that will serve as column names
unique_vars <- unique(unlist(map(strsplit(names(col_means), "subj."), 1)))
unique_IDs <- unique(unlist(map(strsplit(names(col_means), "subj."), 2)))
# Remove putative IDs with "trans" appended in case of z
fake_IDs <- grep("trans", unique_IDs)
if (length(fake_IDs)) unique_IDs <- unique_IDs[-c(fake_IDs)]
map_est <- unlist(foreach(i = 1:ncol(s_df)) %do% {
as.numeric(bayestestR::map_estimate(s_df[, i]))
})
mdf <-
setNames(data.frame(
# Assign the ID as first column..
unique_IDs,
# ..and split into columns for each new variable
matrix(map_est, ncol=length(unique_vars))
),
# Assign column names via setNames
c("ID", unique_vars))
}
mdf$ID <- as.numeric(as.character(mdf$ID))
z_IDs <- grep("trans", mdf$ID)
if (length(z_IDs)) {
mdf[grepl("trans", mdf$ID), "ID"] <- unique(mdf$ID[!grepl("trans", mdf$ID)])
}
# Spot checks that estimates are v similar
#mdf[550:560]
# col_means[550:560]
# mdf$map_est[180:190]
# col_means[180:190]
# mdf$map_est[1:10]
# col_means[1:10]
mdf
}Gelman Rubin of winning model
These files are too big to upload, but the MAP estimates derived from them are in the Public_Data dir — see below
ReadTrace <- function(unique_str) read.csv(unique_str)
ndgr1 <- ReadTrace("../../model_res/final_traces_and_dics/traces/GR_run_also8k_ddm_add_trialwise_NO-SZ_s_vt_poutlier052177.csv")
ndgr2 <- ReadTrace("../../model_res/final_traces_and_dics/traces/GR_run_also8k_ddm_add_trialwise_NO-SZ_s_vt_poutlier052923.csv")
ndgr3 <- ReadTrace("../../model_res/final_traces_and_dics/traces/GR_run_also8k_ddm_add_trialwise_NO-SZ_s_vt_poutlier056294.csv")
ndgr4 <- ReadTrace("../../model_res/final_traces_and_dics/traces/GR_run_also8k_ddm_add_trialwise_NO-SZ_s_vt_poutlier058876.csv")d1c <- rbind(ndgr1, ndgr2, ndgr3, ndgr4)Fine
rhats <- unlist(foreach (i = 2:ncol(ndgr1)) %do% rstan::Rhat(cbind(
unlist(ndgr1[i]),
unlist(ndgr2[i]),
unlist(ndgr3[i]),
unlist(ndgr4[i]))))
hist(rhats)
rhats[which(rhats > 1.1)]## numeric(0)max(rhats)## [1] 1.021465Parameter recovery
# recovered <-
# read.csv("../../model_res/par_recov/HDDM_par_recov_traces_v-val_ddm_add_trialwise_NO-SZ_s_vt_with-z_1239.csv")
recovered <-
read.csv("./Public_Data/HDDM_par_recov_traces_v-val_ddm_add_trialwise_NO-SZ_s_vt_with-z_1239.csv")
#simmed <- read.csv("../../model_res/par_recov/HDDM_sims_for_pr_add_trialwise_NO-SZ_s_vt_poutlier059513.csv")
simmed <- read.csv("./Public_Data/HDDM_sims_for_pr_add_trialwise_NO-SZ_s_vt_poutlier059513.csv")The across-trial variability parameters were run at the group level, so for these we’ll look at the range of means used to generate alongside the recovered posterior
generative_sv <- data.frame(simmed %>% group_by(subj_idx) %>% summarize(msv=mean(sim_sv)))
sv_plot <- ggplot(recovered, aes(x=sv)) +
geom_point(data=generative_sv, aes(x=msv, y=.2), size=3, pch=21, fill="gray57", color="black", alpha=.3) +
geom_point(data=generative_sv, aes(x=mean(msv), y=.2), size=5, pch=21, fill="gray57", color="black", alpha=.7) +
geom_density(color="gray57", fill="white", size=3, alpha=0) +
ga + lp +
xlab("") + ylab("") +
labs(title="Across-trial drift rate") +
#+ #ylab("Higher values = more pos. > neg.") +
theme(axis.text.x = element_text(size=25), axis.ticks=element_blank(), axis.text.y=element_blank()) +
theme(plot.title = element_text(margin=margin(b=-30),
size=12, hjust=.05, vjust=0.1))
sv_plot
# Add st
generative_st <- data.frame(simmed %>% group_by(subj_idx) %>% summarize(mst=mean(sim_st)))
st_plot <- ggplot(recovered, aes(x=st)) +
geom_point(data=generative_st, aes(x=mst, y=.2), size=3, pch=21, fill="gray57", color="black", alpha=.3) +
geom_point(data=generative_st, aes(x=mean(mst), y=.2), size=5, pch=21, fill="gray57", color="black", alpha=.7) +
geom_density(color="gray57", fill="white", size=3, alpha=0) +
ga + lp +
xlab("") + ylab("") +
labs(title="Across-trial non-decision time") +
#+ #ylab("Higher values = more pos. > neg.") +
theme(axis.text.x = element_text(size=25), axis.ticks=element_blank(), axis.text.y=element_blank()) +
theme(plot.title = element_text(margin=margin(b=-30),
size=12, hjust=.05, vjust=0.1))
st_plot
For the rest of the parameters we had subject-level estimates, and we examine the MAP estimates of the recovered traces against the generative subject-level parameters
PlotParRecov <- function(sim_df, opt_df, sim_name, rec_name, parameter_label=NULL) {
sim_df$simulated <- unlist(gen_vals[sim_name])
sim_df$recovered <- unlist(opt_df[rec_name])
cor_res <- cor.test(sim_df$simulated, sim_df$recovered)
lower_lim <- min(sim_df$simulated, sim_df$recovered)
upper_lim <- max(sim_df$simulated, sim_df$recovered)
pr_p <-
ggplot(sim_df, aes(x=simulated, y=recovered)) +
geom_line(aes(x=simulated, y=simulated), size=2) +
geom_point(size=7, alpha=.7, pch=21, fill="gray57") +
ga + ap +
labs(title=paste("r =", round(cor_res$estimate, 2)),
subtitle=parameter_label) +
theme(plot.subtitle = element_text(size = 25)) +
theme(axis.text = element_text(size = 15)) +
tp +
xlim(lower_lim, upper_lim) + ylim(lower_lim, upper_lim)
pr_p
}gen_vals <- simmed %>% group_by(subj_idx) %>% summarize_at(
names(simmed)[grep("sim", names(simmed))], mean
)
rec_map <- GetMapEsts(GetSubjTraces(recovered), "short")
rec_map$conv_z <- InvLogit(rec_map$z_) # Transform recov z back on generative scaleests <- c()
ests[1] <- cor.test(gen_vals$sim_a, rec_map$a_)$est
ests[2] <- cor.test(gen_vals$sim_t, rec_map$t_)$est
ests[3] <- cor.test(gen_vals$sim_z, rec_map$conv_z)$est
ests[4] <- cor.test(gen_vals$sim_vi, rec_map$v_Intercept_)$est
ests[5] <- cor.test(gen_vals$sim_vv, rec_map$v_val_ctr_)$est
ests[6] <- cor.test(gen_vals$sim_vs, rec_map$v_sess_ctr_)$est
ests[7] <- cor.test(gen_vals$sim_vsv, rec_map$v_val_ctr.sess_ctr_)$est
ests## [1] 0.9665570 0.9608423 0.8139605 0.9264843 0.9706376 0.8341225 0.9384575median(ests)## [1] 0.9384575range(ests)## [1] 0.8139605 0.9706376pr_plots <- list()
if (all(rec_map$ID==gen_vals$subj_idx)) {
pr_plots[[1]] <- PlotParRecov(gen_vals, rec_map, "sim_a", "a_", "Threshold (a)")
pr_plots[[2]] <- PlotParRecov(gen_vals, rec_map, "sim_t", "t_", "Non-decision time (t)")
pr_plots[[3]] <- PlotParRecov(gen_vals, rec_map, "sim_z", "conv_z", "Starting point bias (z)")
pr_plots[[4]] <- PlotParRecov(gen_vals, rec_map, "sim_vi", "v_Intercept_", "Drift rate intercept")
pr_plots[[5]] <- PlotParRecov(gen_vals, rec_map, "sim_vv", "v_val_ctr_", "Drift rate valence")
pr_plots[[6]] <- PlotParRecov(gen_vals, rec_map, "sim_vs", "v_sess_ctr_", "Drift rate time")
pr_plots[[7]] <- PlotParRecov(gen_vals, rec_map, "sim_vsv", "v_val_ctr.sess_ctr_", "Drift rate valence * time")
}pr_grid <- (pr_plots[[1]] + pr_plots[[2]]) /
(pr_plots[[3]] + pr_plots[[4]]) /
(pr_plots[[5]] + pr_plots[[6]])
odd_pr_out <- pr_plots[[7]]
pr_grid
#ggsave("../paper/figs/pr-grid_SUPP.png", pr_grid, width=15, height=11, dpi=800)odd_pr_out
#ggsave("../paper/figs/odd-pr_SUPP.png", odd_pr_out, width=7.5, height=11/3, dpi=800)across_trial <- sv_plot + st_plot + plot_annotation(
"Recovered group posterior and generative parameters (mean and individual)",
theme = theme(plot.title = element_text(size = 15)))
across_trial
#ggsave("../paper/figs/across-trial_pr-SUPP.png", across_trial, width=8, height=4)Behavioral results and reflection in modeling
RTs and Endorsements
- Endorsements and RT shift by condition
bx_df$session <- factor(bx_df$session, levels=c("Pre", "Post"))
pre_post_vends <- bx_df %>%
group_by(valence, session) %>%
summarize(mr=mean(response))## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.pre_post_vends_id <- bx_df %>%
group_by(valence, session, subj_idx) %>%
summarize(mr=mean(response))## `summarise()` has grouped output by 'valence', 'session'. You can override
## using the `.groups` argument.pre_post_vends_id## # A tibble: 384 × 4
## # Groups: valence, session [4]
## valence session subj_idx mr
## <chr> <fct> <int> <dbl>
## 1 negative Pre 1 0.556
## 2 negative Pre 2 0.393
## 3 negative Pre 3 0.233
## 4 negative Pre 4 0.103
## 5 negative Pre 6 0.148
## 6 negative Pre 7 0.407
## 7 negative Pre 8 0.423
## 8 negative Pre 9 0.467
## 9 negative Pre 10 0.286
## 10 negative Pre 11 0
## # … with 374 more rowspp_sds <- pre_post_vends_id %>% group_by(valence, session) %>%
summarize(ci95=qnorm(.975)*sd(mr)/sqrt(length(unique(bx_df$subj_idx))))## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.if (all(pp_sds[1:2] == pre_post_vends[1:2])) pre_post_vends$ci95 <- pp_sds$ci95
# Sanity check
#length(unique(bx_df$subj_idx))
# pre_post_vends_id %>% filter(valence == "negative" & session == "Pre") %>%
# select(mr) %>% summarize(sd(mr))/sqrt(96)*qnorm(.975)
#qnorm(.975)#pre_post_vends_id <- bx_df %>% group_by(valence, session, subj_idx) %>% summarize(mr=mean(response))
ends_plot <- ggplot(pre_post_vends_id, aes(x=session, y=mr, group=session, fill=valence)) +
geom_jitter(pch=21, color="black", width=.16, size=3, alpha=.2) +
geom_line(data=pre_post_vends, aes(y=mr, group=valence), color="gray57", pch=21, size=1.5) +
geom_errorbar(data=pre_post_vends, aes(x=session,
ymin=mr-ci95, ymax=mr+ci95), color="black", inherit.aes=FALSE, size=1.5, width=.25) +
geom_point(data=pre_post_vends, aes(y=mr, fill=valence), pch=21, size=6, alpha=.8) +
ylab("proportion endorsements \n (with 95% CI)") + xlab("") +
scale_fill_manual(values=c("red", "blue")) +
ga + ap + tol ## Warning: Ignoring unknown parameters: shapeends_plot
#ggsave("../paper/figs/ends_plot.png", ends_plot, width=11, height=7)pre_post_rts <- bx_df %>% group_by(valence, session) %>% summarize(mr=mean(rt))## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.pre_post_rts_id <- bx_df %>% group_by(valence, session, subj_idx) %>% summarize(mr=mean(rt))## `summarise()` has grouped output by 'valence', 'session'. You can override
## using the `.groups` argument.pprt_sds <- pre_post_rts_id %>% group_by(valence, session) %>%
summarize(ci95=qnorm(.975)*sd(mr)/sqrt(length(unique(bx_df$subj_idx))))## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.if (all(pprt_sds[1:2] == pre_post_rts[1:2])) pre_post_rts$ci95 <- pprt_sds$ci95
rts_plot <- ggplot(pre_post_rts_id, aes(x=session, y=mr, group=session, fill=valence)) +
geom_jitter(pch=21, color="black", width=.16, size=3, alpha=.2) +
geom_line(data=pre_post_rts, aes(y=mr, group=valence), color="gray57", pch=21, size=1.5) +
geom_errorbar(data=pre_post_rts, aes(x=session,
ymin=mr-ci95, ymax=mr+ci95), color="black", inherit.aes=FALSE, size=1.5, width=.25) +
geom_point(data=pre_post_rts, aes(y=mr, fill=valence), pch=21, size=6, alpha=.8) +
ylab("mean reaction time \n (with 95% CI)") + xlab("") +
scale_fill_manual(values=c("red", "blue")) +
ga + ap + tol ## Warning: Ignoring unknown parameters: shaperts_plot
#ggsave("../paper/figs/rt_plot.png", rts_plot, width=11, height=7)Regression models
#scale(bx_df$valence_n)
summary (m1 <- glmer(response ~ val_ctr*sess_ctr +
(1|subj_idx), data=bx_df, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: response ~ val_ctr * sess_ctr + (1 | subj_idx)
## Data: bx_df
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 8184.8 8221.4 -4087.4 8174.8 11051
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4241 -0.3809 0.2075 0.3592 5.1879
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_idx (Intercept) 0.2043 0.452
## Number of obs: 11056, groups: subj_idx, 96
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.08590 0.05474 1.569 0.117
## val_ctr 2.03324 0.03127 65.025 < 2e-16 ***
## sess_ctr -0.02974 0.02943 -1.011 0.312
## val_ctr:sess_ctr 0.17231 0.02944 5.852 4.85e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) vl_ctr sss_ct
## val_ctr 0.038
## sess_ctr 0.065 -0.008
## vl_ctr:sss_ -0.004 0.120 0.072car::vif(m1)## val_ctr sess_ctr val_ctr:sess_ctr
## 1.014929 1.005568 1.020222# Singular - unsuprisingly bc there's just 2 levels of each per subj
# summary (m2 <- glmer(response ~ val_ctr*sess_ctr +
# (val_ctr + sess_ctr|subj_idx), data=bx_df, family="binomial", control = glmerControl(optimizer = "bobyqa")))Using this one bc includes sensible and important RE and not singular
summary (m3 <- glmer(response ~ val_ctr*sess_ctr +
(val_ctr|subj_idx), data=bx_df, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: response ~ val_ctr * sess_ctr + (val_ctr | subj_idx)
## Data: bx_df
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7885.0 7936.2 -3935.5 7871.0 11049
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9890 -0.3662 0.1725 0.3344 4.9228
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subj_idx (Intercept) 0.1699 0.4121
## val_ctr 0.4699 0.6855 -0.08
## Number of obs: 11056, groups: subj_idx, 96
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.07915 0.05488 1.442 0.149
## val_ctr 2.18484 0.07837 27.878 < 2e-16 ***
## sess_ctr -0.03211 0.02975 -1.080 0.280
## val_ctr:sess_ctr 0.18136 0.02976 6.094 1.1e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) vl_ctr sss_ct
## val_ctr -0.039
## sess_ctr 0.073 -0.004
## vl_ctr:sss_ -0.005 0.051 0.069car::vif(m3)## val_ctr sess_ctr val_ctr:sess_ctr
## 1.00267 1.00484 1.00746summary (rm1 <- lmer(rt ~ val_ctr*sess_ctr +
(1|subj_idx), data=bx_df))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ val_ctr * sess_ctr + (1 | subj_idx)
## Data: bx_df
##
## REML criterion at convergence: 12762.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2574 -0.7110 -0.2027 0.5462 4.4470
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_idx (Intercept) 0.04649 0.2156
## Residual 0.17980 0.4240
## Number of obs: 11056, groups: subj_idx, 96
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.520e+00 2.237e-02 9.485e+01 67.954 <2e-16 ***
## val_ctr -3.703e-02 4.034e-03 1.096e+04 -9.180 <2e-16 ***
## sess_ctr -3.649e-02 4.033e-03 1.096e+04 -9.048 <2e-16 ***
## val_ctr:sess_ctr 5.267e-03 4.034e-03 1.096e+04 1.306 0.192
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) vl_ctr sss_ct
## val_ctr 0.001
## sess_ctr 0.000 0.001
## vl_ctr:sss_ 0.000 -0.001 0.006car::vif(rm1)## val_ctr sess_ctr val_ctr:sess_ctr
## 1.000001 1.000031 1.000031Used for same reason
summary (rm2 <- lmer(rt ~ val_ctr*sess_ctr +
(val_ctr|subj_idx), data=bx_df))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ val_ctr * sess_ctr + (val_ctr | subj_idx)
## Data: bx_df
##
## REML criterion at convergence: 12701.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3659 -0.7089 -0.2026 0.5390 4.5891
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subj_idx (Intercept) 0.046590 0.21585
## val_ctr 0.002408 0.04907 0.24
## Residual 0.177400 0.42119
## Number of obs: 11056, groups: subj_idx, 96
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.520e+00 2.239e-02 9.485e+01 67.901 < 2e-16 ***
## val_ctr -3.690e-02 6.414e-03 9.474e+01 -5.752 1.08e-07 ***
## sess_ctr -3.640e-02 4.007e-03 1.086e+04 -9.084 < 2e-16 ***
## val_ctr:sess_ctr 5.124e-03 4.007e-03 1.086e+04 1.279 0.201
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) vl_ctr sss_ct
## val_ctr 0.188
## sess_ctr 0.000 0.000
## vl_ctr:sss_ 0.000 -0.001 0.006car::vif(rm1)## val_ctr sess_ctr val_ctr:sess_ctr
## 1.000001 1.000031 1.000031DIC
model_names <- c("baseline", "+ valence", "+ valence * time", "+ bias", "+ sv and st")
dic_values <-
c(
# read.csv("./../../model_res/final_traces_and_dics/dic/s_b_389_dic.csv", header=FALSE)$V2,
# read.csv("./../../model_res/final_traces_and_dics/dic/s_b1_3987_dic.csv", header=FALSE)$V2,
# read.csv("./../../model_res/final_traces_and_dics/dic/s_vt_dic.csv", header=FALSE)$V2,
# read.csv("./../../model_res/final_traces_and_dics/dic/GR_8k_s_vt_poutlier058886dic.csv", header=FALSE)$V2,
# read.csv("./../../model_res/final_traces_and_dics/dic/GR_run_also8k_ddm_add_trialwise_NO-SZ_s_vt_poutlier056294dic.csv", header=FALSE)$V2)
read.csv("./Public_Data/dic/s_b_389_dic.csv", header=FALSE)$V2,
read.csv("./Public_Data/dic/s_b1_3987_dic.csv", header=FALSE)$V2,
read.csv("./Public_Data/dic/s_vt_dic.csv", header=FALSE)$V2,
read.csv("./Public_Data/dic/GR_8k_s_vt_poutlier058886dic.csv", header=FALSE)$V2,
read.csv("./Public_Data/dic/GR_run_also8k_ddm_add_trialwise_NO-SZ_s_vt_poutlier056294dic.csv", header=FALSE)$V2)
dic_df <- data.frame(model_names, dic_values)
dic_df$delta_dic <- c(0, dic_values[2:5]-dic_values[1])
#dic_df$delta_dic <- as.numeric(c(0, unlist(dic_df[2:nrow(dic_df), 1]) - unlist(dic_df[1, 1])))
dic_df$model_names <- factor(dic_df$model_names, levels=c("baseline", "+ valence", "+ valence * time", "+ bias", "+ sv and st"))
dic_p <- ggplot(dic_df[2:5, ], aes(x=model_names, y=delta_dic)) +
geom_bar(stat="identity", color="black", fill="white", size=2) +
ga + ap + tol +
ylab(TeX("$\\Delta$ DIC")) +
xlab("") +
labs(title = "Model comparison",
subtitle = "Relative to baseline (no regressors) model") +
theme(plot.subtitle = element_text(size = 15)) +
#ggtitle("Model comparison relative to baseline") +
tp +
theme(axis.text.x = element_text(angle=30, hjust=1, size=20))#+ xlim(500, -1200)
dic_p
#ggsave("../paper/figs/dic.png", dic_p, width=8, height=6, dpi=600)Posteriors
val_posteriors <- ggplot(d1c, aes(x=v_val_ctr)) +
geom_vline(xintercept=0, size=2) +
geom_density(color="chocolate", fill="white", alpha=1, size=3) +
geom_density(aes(x=v_val_ctr.sess_ctr), color="tan1", fill="white", alpha=1, size=3) +
xlim(-.1, 2) + #geom_vline(xintercept = 0, color="gray57", size=2)
ga + lp +
xlab("") + ylab("") + #ylab("Higher values = more pos. > neg.") +
labs(title="Drift valence (dark) and valence*time (light) \n regressor posteriors") +
theme(axis.text.x = element_text(size=35), axis.ticks=element_blank(), axis.text.y=element_blank()) +
theme(plot.title = element_text(margin=margin(b=-30),
size=30, hjust=.18, vjust=.6))
val_posteriors
#ggsave("../paper/figs/val_posteriors.png", val_posteriors, width=12, height=7)bayestestR::map_estimate(d1c$v_val_ctr)## MAP Estimate: 1.33which(d1c$v_val_ctr < 0)## integer(0)#hist(d1c$v_val_ctr)
bayestestR::map_estimate(d1c$v_val_ctr.sess_ctr)## MAP Estimate: 0.13which(d1c$v_val_ctr.sess_ctr < 0)## integer(0)#hist(d1c$v_val_ctr.sess_ctr)Endorsement proportions
Positive and negative pre and differences therein
nends_b <- bx_df %>% filter(session=="Pre", valence=="negative") %>%
group_by(subj_idx) %>% summarize(m=mean(response), mrt=mean(rt))
pends_b <- bx_df %>% filter(session=="Pre", valence=="positive") %>% group_by(subj_idx) %>%
summarize(m=mean(response), mrt=mean(rt))
if (all(pends_b$subj_idx==nends_b$subj_idx)) {
pends_b$pn_diff <- pends_b$m-nends_b$m
pends_b$pn_rt_diff <- pends_b$mrt-nends_b$mrt
}
# pends_b$mrt
# nends_b$mrtRTs
rt <- bx_df %>%
group_by(subj_idx) %>% summarize(m=mean(rt), mrt=mean(rt))
rt_val <- bx_df %>%
group_by(subj_idx, valence) %>% summarize(mrt=mean(rt))## `summarise()` has grouped output by 'subj_idx'. You can override using the
## `.groups` argument.rt_val_short <- data.frame(
"ID"=unique(rt_val$subj_idx),
"overall_rt_diff"=rt_val[rt_val$valence == "positive", "mrt"]-rt_val[rt_val$valence == "negative", "mrt"])Endorsements at post and differences from pre to post
nends_p <- bx_df %>% filter(session=="Post", valence=="negative") %>%
group_by(subj_idx) %>% summarize(m=mean(response), mrt=mean(rt))
pends_p <- bx_df %>% filter(session=="Post", valence=="positive") %>% group_by(subj_idx) %>% summarize(m=mean(response), mrt=mean(rt))
if (all(pends_p$subj_idx==nends_p$subj_idx)) {
pends_p$pn_diff <- pends_p$m-nends_p$m
pends_p$pn_rt_diff <- pends_p$mrt-nends_p$mrt
}d1m <- GetMapEsts(GetSubjTraces(d1c), "short")The traces are massive so can’t put these in the public data folder, however have uploaded the map ests — thus these can be read in to reproduce further analysis
#write.csv(d1m, "./Public_Data/short_form_map_ests.csv")
d1m <- read.csv("./Public_Data/short_form_map_ests.csv")Write
rt <- bx_df %>%
group_by(subj_idx) %>% summarize(m=mean(rt), mrt=mean(rt))
if (all(d1m$ID==pends_p$subj_idx) & all(d1m$ID == pends_b$subj_idx)) {
cat("\n DDM baseline ends"); print(cor.test(d1m$v_val_ctr_, pends_b$pn_diff))
cat("\n DDM baseline rt"); print(cor.test(d1m$v_val_ctr_, pends_b$pn_rt_diff))
cat("\n DDM change ends");print(cor.test(d1m$v_val_ctr.sess_ctr_, pends_p$pn_diff-pends_b$pn_diff))
cat("\n DDM change rt"); print(cor.test(d1m$v_val_ctr.sess_ctr_, pends_p$pn_rt_diff-pends_b$pn_rt_diff))
} ##
## DDM baseline ends
## Pearson's product-moment correlation
##
## data: d1m$v_val_ctr_ and pends_b$pn_diff
## t = 12.389, df = 94, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6971071 0.8532895
## sample estimates:
## cor
## 0.787528
##
##
## DDM baseline rt
## Pearson's product-moment correlation
##
## data: d1m$v_val_ctr_ and pends_b$pn_rt_diff
## t = -0.61816, df = 94, p-value = 0.538
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2607878 0.1386257
## sample estimates:
## cor
## -0.06362865
##
##
## DDM change ends
## Pearson's product-moment correlation
##
## data: d1m$v_val_ctr.sess_ctr_ and pends_p$pn_diff - pends_b$pn_diff
## t = 8.11, df = 94, p-value = 1.876e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5062421 0.7461185
## sample estimates:
## cor
## 0.6416083
##
##
## DDM change rt
## Pearson's product-moment correlation
##
## data: d1m$v_val_ctr.sess_ctr_ and pends_p$pn_rt_diff - pends_b$pn_rt_diff
## t = 0.61594, df = 94, p-value = 0.5394
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1388492 0.2605754
## sample estimates:
## cor
## 0.06340168# This is the overall RT valence RT averaging over time point
if (all(d1m$ID == rt_val_short$ID)) cat("\n DDM val rt"); print(cor.test(d1m$v_val_ctr_, rt_val_short$mrt))##
## DDM val rt##
## Pearson's product-moment correlation
##
## data: d1m$v_val_ctr_ and rt_val_short$mrt
## t = 0.53113, df = 94, p-value = 0.5966
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1474027 0.2524175
## sample estimates:
## cor
## 0.05469969if (all(d1m$ID == rt$subj_idx)) cat("\n DDM threshold - rt"); print(cor.test(d1m$a_, rt$m))##
## DDM threshold - rt##
## Pearson's product-moment correlation
##
## data: d1m$a_ and rt$m
## t = 9.2603, df = 94, p-value = 6.877e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5690004 0.7827845
## sample estimates:
## cor
## 0.6906943Model validation w/ PPCs
s_bdf <- bx_df#m1_ddm_ppcs <- read.csv("../../model_res/ppcs/HDDM_ppcs_v-val_ddm_add_trialwise_NO-SZ_s_vt_with-z_9533.csv")
m1_ddm_ppcs <- read.csv("./Public_Data/HDDM_ppcs_v-val_ddm_add_trialwise_NO-SZ_s_vt_with-z_9533.csv")Recode response
m1_ddm_ppcs[m1_ddm_ppcs$sess_ctr < 0, "session"] <- "Pre"
m1_ddm_ppcs[m1_ddm_ppcs$sess_ctr > 0, "session"] <- "Post"
m1_ddm_ppcs$session <- factor(m1_ddm_ppcs$session, levels=c("Pre", "Post"))
m1_ddm_ppcs[m1_ddm_ppcs$val_ctr > 0, "valence"] <- "positive"
m1_ddm_ppcs[m1_ddm_ppcs$val_ctr < 0, "valence"] <- "negative"
m1_ddm_ppcs$valence <- factor(m1_ddm_ppcs$valence, levels=c("negative", "positive"))
#qp_no_t$response <- factor(qp_no_t$response, levels=c("no", "yes"))hist(m1_ddm_ppcs$rt, breaks = 100, xlim = c(0, 4))
Quantiles look pretty good
quantile(s_bdf$rt, probs = seq(0.05, .95, .1))## 5% 15% 25% 35% 45% 55% 65% 75%
## 0.901915 1.045560 1.151572 1.255880 1.360530 1.483183 1.627107 1.799898
## 85% 95%
## 2.042533 2.460470quantile(m1_ddm_ppcs$rt, probs = seq(0.05, .95, .1))## 5% 15% 25% 35% 45% 55% 65% 75%
## 0.8895856 1.0505985 1.1580293 1.2532260 1.3480761 1.4507999 1.5748151 1.7384998
## 85% 95%
## 1.9930176 2.5840549Code flip rt for PPC plots
s_bdf$flip_rt <- s_bdf$rt
s_bdf[s_bdf$response==0, "flip_rt"] <- -s_bdf[s_bdf$response==0, "flip_rt"]
m1_ddm_ppcs$flip_rt <- m1_ddm_ppcs$rt
m1_ddm_ppcs[m1_ddm_ppcs$response==0, "flip_rt"] <- -m1_ddm_ppcs[m1_ddm_ppcs$response==0, "flip_rt"]
m1_ddm_ppcs$flip_rt <- m1_ddm_ppcs$rt
m1_ddm_ppcs[m1_ddm_ppcs$response==0, "flip_rt"] <- -m1_ddm_ppcs[m1_ddm_ppcs$response==0, "flip_rt"] Cap the PPCs at 3 because that’s where the task cut off
m1_ddm_ppcs <- m1_ddm_ppcs %>% filter(rt < 3)Create alt negative coded rt with negative side flipped so lines up with QP plots
m1_ddm_ppcs$alt_rt <- m1_ddm_ppcs$flip_rt
s_bdf$alt_rt <- s_bdf$flip_rt
#m1_ddm_ppcs
m1_ddm_ppcs[m1_ddm_ppcs$valence=="negative", "alt_rt"] <- -m1_ddm_ppcs[m1_ddm_ppcs$valence=="negative", "alt_rt"]
s_bdf[s_bdf$valence=="negative", "alt_rt"] <- -s_bdf[s_bdf$valence=="negative", "alt_rt"]ddm_alt <- ggplot(s_bdf, aes(x=alt_rt)) +
geom_density(alpha=.1, size=1.3, color="black", fill="black") +
geom_density(data=m1_ddm_ppcs, alpha=0, color="darkorange",
size=1.3, fill="red", linetype="longdash") +
facet_wrap( ~ valence) +
ga + ap + ft + tp + xlab("reaction time") + xlim(-3.1, 3.1) +
theme(axis.text=element_blank(), axis.ticks = element_blank()) +
ylab("")
ddm_alt
#ggsave("../paper/figs/big_ppc_alt.png", ddm_alt, width=9.25, height=3)X axis - response frequency
ddm_neg <-
table(m1_ddm_ppcs[m1_ddm_ppcs$valence=="negative", "response"])[2]/
sum(table(m1_ddm_ppcs[m1_ddm_ppcs$valence=="negative", "response"]))
ddm_pos <-
table(m1_ddm_ppcs[m1_ddm_ppcs$valence=="positive", "response"])[2]/
sum(table(m1_ddm_ppcs[m1_ddm_ppcs$valence=="positive", "response"]))
emp_neg <-
table(s_bdf[s_bdf$valence=="negative", "response"])[2]/
sum(table(s_bdf[s_bdf$valence=="negative", "response"]))
emp_pos <-
table(s_bdf[s_bdf$valence=="positive", "response"])[2]/
sum(table(s_bdf[s_bdf$valence=="positive", "response"]))Y axis - reaction time quantiles
qs <- seq(0.1, .9, .2)CalculateRTQuantile <- function(val_resp_df, qs) {
### Get subject average quantiles at a given valence-response level ###
all_subj_quantiles <-
lapply(split(val_resp_df, val_resp_df$subj_idx), function(x) {
out <-
data.table(unique(x$subj_idx),
unique(x$valence),
unique(x$response),
quantile(x$rt, probs = qs),
qs) %>%
setNames(c("ID", "valence", "response", "quantiles", "qs"))
out
}) %>% bind_rows()
all_subj_quantiles
}ddm <-
data.table(
rbind(
CalculateRTQuantile(m1_ddm_ppcs %>% filter(valence=="negative", response==1), qs),
CalculateRTQuantile(m1_ddm_ppcs %>% filter(valence=="positive", response==1), qs),
CalculateRTQuantile(m1_ddm_ppcs %>% filter(valence=="negative", response==0), qs),
CalculateRTQuantile(m1_ddm_ppcs %>% filter(valence=="positive", response==0), qs)),
"model"="ddm"
)
emp <- data.table(
rbind(
CalculateRTQuantile(s_bdf %>% filter(valence=="negative", response==1), qs),
CalculateRTQuantile(s_bdf %>% filter(valence=="positive", response==1), qs),
CalculateRTQuantile(s_bdf %>% filter(valence=="negative", response==0), qs),
CalculateRTQuantile(s_bdf %>% filter(valence=="positive", response==0), qs)),
"model"="empirical")
all_subj_qs <- rbind(ddm, emp)Put it all together
q_rt_sum <- all_subj_qs %>% group_by(model, valence, response, qs) %>% summarize(q=mean(quantiles))## `summarise()` has grouped output by 'model', 'valence', 'response'. You can
## override using the `.groups` argument.#q_rt_sum %>% filter(response==1 && valence=="negative", model=="ddm")
#ddm_neg
qp_df <-
rbind(
## DDM
# For each model-valence, use the RT quantiles at response 1 and 0, whose response probs are p(endorse) and 1-p(endorse)
q_rt_sum %>% filter(response==1 && valence=="negative", model=="ddm") %>% mutate(rp=ddm_neg),
q_rt_sum %>% filter(response==0 && valence=="negative", model=="ddm") %>% mutate(rp=1-ddm_neg),
q_rt_sum %>% filter(response==1 && valence=="positive", model=="ddm") %>% mutate(rp=ddm_pos),
q_rt_sum %>% filter(response==0 && valence=="positive", model=="ddm") %>% mutate(rp=1-ddm_pos),
## Empirical
q_rt_sum %>% filter(response==1 && valence=="negative", model=="empirical") %>% mutate(rp=emp_neg),
q_rt_sum %>% filter(response==0 && valence=="negative", model=="empirical") %>% mutate(rp=1-emp_neg),
q_rt_sum %>% filter(response==1 && valence=="positive", model=="empirical") %>% mutate(rp=emp_pos),
q_rt_sum %>% filter(response==0 && valence=="positive", model=="empirical") %>% mutate(rp=1-emp_pos)
)qp_df$model <- factor(qp_df$model, levels=c("empirical", "ddm"))qp_plot <- ggplot(qp_df, aes(x=rp, y=q, fill=model)) + geom_point(pch=21, size=6, alpha=.7) +
facet_wrap( ~ valence) + ga + ap + lp + ft +
ylab("reaction time") + xlab("response frequency") +
scale_fill_manual(values=c("black", "darkorange"), labels=c("Empirical", "DDM")) +
theme(axis.text.x = element_text(angle=30, hjust=1, size=20)) +
theme(legend.position = c(.25, .65)) #+ xlim(500, -1200)
qp_plot
#ggsave("../paper/figs/qp_plot.png", qp_plot, width=10, height=4.25)# probably don't need this plot
# ppc_lf <- rbind(
# data.table(m1_ddm_ppcs %>% select("subj_idx", "response", "valence", "flip_rt"), "ddm"),
# data.table(lm1_ppc %>% select("subj_idx", "response"="responses", "valence", "flip_rt"), "levy"),
# data.table(s_bdf %>% select("subj_idx", "response", "valence", "flip_rt"), "empirical")
# ) %>% setNames(c("ID", "response", "valence", "flip_rt", "model"))
#
# table(ppc_lf$response)
#
# # Summarize neg yes
# ny <- ppc_lf %>% filter( valence=="negative") %>% group_by(ID, model) %>% summarize(mr=mean(response))
#
#
# sd(as.numeric(unlist(ny[ny$model=="ddm", "mr"])))/sqrt(length(unique(bx_df$subj_idx)))
# nys <- ny %>% group_by(model) %>% summarize(m=mean(mr))
#
# nys$se <-
# c(
# sd(as.numeric(unlist(ny[ny$model=="ddm", "mr"])))/sqrt(length(unique(bx_df$subj_idx))),
# sd(as.numeric(unlist(ny[ny$model=="levy", "mr"])))/sqrt(length(unique(bx_df$subj_idx))),
# sd(as.numeric(unlist(ny[ny$model=="empirical", "mr"])))/sqrt(length(unique(bx_df$subj_idx)))
# )
#
# # c(
# # sd(unlist(ppc_lf[ppc_lf$model=="ddm", "response"]))/sqrt(length(unique(bx_df$subj_idx))), # ** make sure right
# # sd(unlist(ppc_lf[ppc_lf$model=="levy", "response"]))/sqrt(length(unique(bx_df$subj_idx))),
# # sd(unlist(ppc_lf[ppc_lf$model=="empirical", "response"]))/sqrt(length(unique(bx_df$subj_idx))))
#
#
# pd <- .25
# ny$model <- factor(ny$model, levels=c("ddm", "empirical", "levy"))
# neg <- ggplot(ny, aes(x=model, y=mr, group=ID, fill=model)) +
# #geom_line(position=position_dodge(width=.1), alpha=.16) +
# geom_point(width=.1, position=position_dodge(width=.18), alpha=1, size = 2, pch=21) +
# scale_fill_manual(values=c("orange", "gray57", "purple")) +
# # geom_errorbar(data=nys, aes(x=model, ymax=m+se, ymin=m-se), inherit.aes = FALSE, size=2.2, width=.15, color="blue") +
# geom_point(data=nys, aes(x=model, y=m), inherit.aes = FALSE, size=8, pch=21, fill="white") +
#
# geom_hline(yintercept = c(unlist(nys[nys$model=="empirical", "m"])), size=1.2) +
# ylim(0, 1) + ga + ap + tol + xlab("") + ylab("response \n frequency") +
# theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title = element_text(size=55))
#
# neg
#
# nys
# Summarize neg yes
# ny <- ppc_lf %>% filter( valence=="positive") %>% group_by(ID, model) %>% summarize(mr=mean(response))
#
# nys <- ny %>% group_by(model) %>% summarize(m=mean(mr))
#
# nys$se <-
# c(
# sd(as.numeric(unlist(ny[ny$model=="ddm", "mr"])))/sqrt(length(unique(bx_df$subj_idx))),
# sd(as.numeric(unlist(ny[ny$model=="levy", "mr"])))/sqrt(length(unique(bx_df$subj_idx))),
# sd(as.numeric(unlist(ny[ny$model=="empirical", "mr"])))/sqrt(length(unique(bx_df$subj_idx)))
# )
# #
# # c(
# # sd(unlist(ppc_lf[ppc_lf$model=="ddm", "response"]))/sqrt(length(unique(bx_df$subj_idx))), # ** make sure right
# # sd(unlist(ppc_lf[ppc_lf$model=="levy", "response"]))/sqrt(length(unique(bx_df$subj_idx))),
# # sd(unlist(ppc_lf[ppc_lf$model=="empirical", "response"]))/sqrt(length(unique(bx_df$subj_idx))))
#
#
# pd <- .25
# ny$model <- factor(ny$model, levels=c("ddm", "empirical", "levy"))
# pos <- ggplot(ny, aes(x=model, y=mr, group=ID, fill=model)) +
# geom_point(width=.1, position=position_dodge(width=.18), alpha=1, size = 2, pch=21) +
# scale_fill_manual(values=c("orange", "gray57", "purple")) +
# # geom_errorbar(data=nys, aes(x=model, ymax=m+se, ymin=m-se), inherit.aes = FALSE, size=2.2, width=.15, color="blue") +
# geom_point(data=nys, aes(x=model, y=m), inherit.aes = FALSE, size=8, pch=21, fill="white") +
#
# geom_hline(yintercept = c(unlist(nys[nys$model=="empirical", "m"])), size=1.2) +
# ylim(0, 1) + ga + ap + tol + xlab("") + ylab("") +
# theme(axis.text=element_blank(), axis.ticks=element_blank(), axis.title = element_text(size=40)) #+
# #scale_fill_manual(values=c("black", "orange", "purple"), labels=c("Empirical", "DDM", "Lévy"))
# #+ xlim(500, -1200)
#
# #pos
# both <- neg+pos
# both
#ggsave("../paper/figs/ppc_freqs.png", both, width=18, height=6)Symptom analyses and regression models
First prep questionnaire items
completer_ids <- unique(bx_df$subj_idx)
#item_qairre_df <- read.csv("./../../data/raw_files/justDASSandFFMQ_with_imputes_replaced_with_999.csv")
item_qairre_df <- read.csv("./Public_Data/justDASSandFFMQ_with_imputes_replaced_with_999.csv")
#all(item_qairre_df$ID.ffmq==item_qairre_df$ID.DASS) # Pulled from different excel sheets so this double check IDs lined up exactly# No reverse scoring
dass_dep_qs <- c(3, 5, 10, 13, 16, 17, 21, 24, 26, 31, 34, 37, 38, 42)
iq_c <- item_qairre_df %>% filter(ID.ffmq %in% completer_ids)DASS
iq_c_d <- iq_c[grep("DASS", names(iq_c))]Time 1 - just depression subscale
# Time 1
dd1 <- iq_c_d[paste0("DASS1.", dass_dep_qs)]# 1.26% of baseline data missing
sum(length(which(dd1==999)), length(which(is.na(dd1))))/(dim(dd1)[1]*dim(dd1)[2])*100## [1] 1.264881dd1 <- data.frame(dd1)
dd1[dd1==999] <- NA
which(is.na(rowSums(dd1))) # 4 pts have NAs ## [1] 2 42 79 83# 1 all data missing, the other just a couple for each
dd1[which(is.na(rowSums(dd1))), ]## DASS1.3 DASS1.5 DASS1.10 DASS1.13 DASS1.16 DASS1.17 DASS1.21 DASS1.24
## 2 2 1 3 2 1 2 2 2
## 42 NA NA NA NA NA NA NA NA
## 79 1 1 1 NA 1 0 0 0
## 83 2 3 3 2 2 1 1 2
## DASS1.26 DASS1.31 DASS1.34 DASS1.37 DASS1.38 DASS1.42
## 2 2 1 1 1 1 NA
## 42 NA NA NA NA NA NA
## 79 1 0 0 0 0 0
## 83 2 2 NA 1 1 3# Put in ID
dd1 <- data.frame("ID"=iq_c_d$ID.DASS, dd1)
# For robustness check, casewise delete pts with any data missing
dd1_casewise_delete <- dd1[-which(is.na(rowSums(dd1))), ]
dd1_casewise_delete$dass_sum_score <- rowSums(dd1_casewise_delete[2:15])
#dd1[which((is.na(rowSums(dd1[2:15])))), ]
id_to_del <- 50 # This is the pt with all missing—since have no data, casewise delete this person
dd1 %>% filter(ID == id_to_del)## ID DASS1.3 DASS1.5 DASS1.10 DASS1.13 DASS1.16 DASS1.17 DASS1.21 DASS1.24
## 1 50 NA NA NA NA NA NA NA NA
## DASS1.26 DASS1.31 DASS1.34 DASS1.37 DASS1.38 DASS1.42
## 1 NA NA NA NA NA NAdd1 <- dd1[!dd1$ID==50, ]
# .21% of items missing after excluding all missing
length(which(is.na(dd1)))/(dim(dd1)[1]*dim(dd1)[2])*100## [1] 0.2105263# Replace any NAs with the mean of the remainding items
for (r in 1:nrow(dd1)) {
if (any(is.na(dd1[r, 2:15]))) {
this_row <- as.numeric(unlist(dd1[r, 2:15])) # Exclude the ID
# Mean to replace NAs with
mean_non_na <- mean(this_row[which(!is.na(this_row))])
this_row[is.na(as.numeric(unlist(this_row)))] <- mean_non_na
# Put in df
dd1[r, 2:15] <- this_row
}
}
dd1$dass_sum_score <- rowSums(dd1[2:15])
#gdf <- haven::read_sav("./../data/raw_files/Britton_ConcurrentStudyA_7.11.16.sav")
#gdf$DASS_Depression_SUM_Q1==rowSums(dd1)mean(dd1$dass_sum_score) # Depression only ## [1] 9.578138sd(dd1$dass_sum_score)## [1] 7.694246Time 2
# Time 2
# 2.31 % of final data missing
dd2 <- iq_c_d[paste0("DASS2.", dass_dep_qs)]
sum(length(which(dd2==999)), length(which(is.na(dd2))))/ (dim(dd2)[1]*dim(dd2)[2])*100## [1] 2.306548dd2[dd2==999] <- NA
which(is.na(rowSums(dd2))) # 3pts have NAs ## [1] 48 50 51# 1 all data missing, the other just a couple for each
dd2[which(is.na(rowSums(dd2))), ]## DASS2.3 DASS2.5 DASS2.10 DASS2.13 DASS2.16 DASS2.17 DASS2.21 DASS2.24
## 48 NA NA NA NA NA NA NA NA
## 50 NA NA NA NA NA NA NA NA
## 51 0 NA 1 NA 1 0 1 1
## DASS2.26 DASS2.31 DASS2.34 DASS2.37 DASS2.38 DASS2.42
## 48 NA NA NA NA NA NA
## 50 NA NA NA NA NA NA
## 51 1 1 NA 1 0 0# Put in ID
dd2 <- data.frame("ID"=iq_c_d$ID.DASS, dd2)
# For robustness check, casewise delete pts with any data missing
dd2_casewise_delete <- dd2[-which(is.na(rowSums(dd2))), ]
dd2_casewise_delete$dass_sum_score <- rowSums(dd2_casewise_delete[2:15])
# Delete pts with all data missing
ids_to_del <- c(56, 58)
dd2 <- dd2 %>% filter(!ID %in% ids_to_del)
# .21% of items missing
length(which(is.na(dd2)))/(dim(dd2)[1]*dim(dd2)[2])*100## [1] 0.212766# Replace the one person missing just a couple items with the mean
mean_to_repl <- mean(unlist(dd2[dd2$ID==59, 2:15][which(!is.na(dd2[dd2$ID==59, 2:15]))]))
dd2[dd2$ID==59, which(is.na(dd2[dd2$ID==59, ]))] <- mean_to_repl
dd2$dass_sum_score <- rowSums(dd2[2:15])
dd2$dass_sum_score## [1] 5.000000 17.000000 1.000000 8.000000 1.000000 2.000000 3.000000
## [8] 0.000000 4.000000 13.000000 3.000000 0.000000 24.000000 1.000000
## [15] 2.000000 1.000000 10.000000 0.000000 0.000000 6.000000 0.000000
## [22] 0.000000 3.000000 1.000000 0.000000 0.000000 0.000000 11.000000
## [29] 0.000000 2.000000 8.000000 0.000000 5.000000 8.000000 2.000000
## [36] 1.000000 18.000000 16.000000 2.000000 14.000000 2.000000 2.000000
## [43] 4.000000 11.000000 1.000000 1.000000 5.000000 0.000000 8.909091
## [50] 6.000000 0.000000 1.000000 3.000000 0.000000 4.000000 8.000000
## [57] 2.000000 1.000000 2.000000 0.000000 19.000000 0.000000 6.000000
## [64] 15.000000 1.000000 22.000000 4.000000 5.000000 5.000000 0.000000
## [71] 0.000000 1.000000 2.000000 4.000000 1.000000 4.000000 0.000000
## [78] 8.000000 2.000000 0.000000 7.000000 3.000000 9.000000 14.000000
## [85] 1.000000 8.000000 6.000000 4.000000 1.000000 0.000000 0.000000
## [92] 2.000000 3.000000 2.000000#any(is.na(dd2))Symptom analyses (which just involve depression questionnaire)
dass_time_1 <- dd1
dass_time_2 <- dd2Make sure the DASS are numerically ordered, then give them names
order <- as.numeric(unlist(
lapply(names(dass_time_1)[grep("DASS", names(dass_time_1))], function(x) unlist(map(strsplit(x, "[.]"), 2)))
))
if (all(unlist(foreach(i = 1:(length(order)-1)) %do% { order[i] < order[i+1]}))) {
extra_d1 <- dass_time_1[names(dass_time_1)[grep("DASS", names(dass_time_1))]] %>%
setNames(c(
# 3, 5, 10, 13
"no-positive-feelings", "couldnt-get-going", "easily-upset", "sad-and-depressed",
# 16, 17, 21
"lost-interest", "not-worth-much-as-person", "life-not-worthwhile",
# 24, 26, 31, 34, 37,
"no-enjoyment", "down-hearted-and-blue", "no-enthusiasm", "worthless", "no-hope",
#38, 42
"life-meaningless", "no-initiative"
))
dass_time_1 <- data.frame(extra_d1, dass_time_1)
}
# Spot checks
#dass_time_1$no.initiative==dd1$DASS1.42
# dass_time_1$no.enjoyment==dd1$DASS1.24
# dass_time_1$life.meaningless==dd1$DASS1.38
# dass_time_1$sad.and.depressed==dd1$DASS1.13order2 <- as.numeric(unlist(
lapply(names(dass_time_2)[grep("DASS", names(dass_time_2))], function(x) unlist(map(strsplit(x, "[.]"), 2)))
))
if (all(unlist(foreach(i = 1:(length(order2)-1)) %do% { order[i] < order2[i+1]}))) {
extra_d2 <- dass_time_2[names(dass_time_2)[grep("DASS", names(dass_time_2))]] %>%
setNames(c(
# 3, 5, 10, 13
"no-positive-feelings", "couldnt-get-going", "easily-upset", "sad-and-depressed",
# 16, 17, 21
"lost-interest", "not-worth-much-as-person", "life-not-worthwhile",
# 24, 26, 31, 34, 37,
"no-enjoyment", "down-hearted-and-blue", "no-enthusiasm", "worthless", "no-hope",
#38, 42
"life-meaningless", "no-initiative"
))
dass_time_2 <- data.frame(extra_d2, dass_time_2)
}
# Spot checks
#dass_time_2$no.initiative==dd2$DASS2.42
#dass_time_2$down.hearted.and.blue==dd2$DASS2.26
#dass_time_2$no.positive.feelings==dd2$DASS2.3
#dass_time_2$sad.and.depressed==dd2$DASS2.13Create a long-form version
d1_short <- dass_time_1 %>% filter(ID %in% dass_time_2$ID)
d1_lf <- d1_short[1:(length(dass_dep_qs)+1)] %>%
pivot_longer(cols = no.positive.feelings:no.initiative,
names_to = "item",
values_to = "scores"
)d2_short <- dass_time_2 %>% filter(ID %in% dass_time_1$ID)
if (all(d1_short$ID==d2_short$ID)) {
change_extra <- d2_short[1:length(dass_dep_qs)] - d1_short[1:length(dass_dep_qs)]
change_df <- data.frame(
change_extra, change_extra %>% setNames(paste0("change_", names(change_extra))))
}
# Changes the order so not working yet
labs <- unlist(lapply(strsplit(unique(d1_lf$item), "[.]"), function(x) {
out <- paste0(x, collapse=" ")
out
}))Drop the few imputed ones for vis
d1_lf <- d1_lf[-c(which(d1_lf$scores %% 1 != 0)), ]Baseline
d1_summs <- d1_lf %>% group_by(item) %>% summarize(m=mean(scores))
d1_summs_ID <- d1_lf %>% group_by(item, ID) %>% summarize(m=mean(scores)) ## `summarise()` has grouped output by 'item'. You can override using the
## `.groups` argument.d1_sds <- d1_lf %>% group_by(item) %>%
summarize(ci95=qnorm(.975)*sd(scores)/sqrt(length(unique(bx_df$subj_idx))))
if (all(d1_summs$item==d1_sds$item)) d1_summs$ci95 <- d1_sds$ci95
d1_summs_arr <-
d1_summs %>% arrange(m)
#d1_summs_arr$item_f <- factor(ch_summs_arr$item, levels=rev(ch_summs_arr$item))
d1_lf$item_f <- factor(d1_lf$item, levels=c(d1_summs_arr$item))
base_labs <- as.character(levels(d1_lf$item_f))
#unlist(map(strsplit(levels(d1_lf$item_f), "change_"), 2))
b_labs <- unlist(lapply(strsplit(base_labs, "[.]"), function(x) {
out <- paste0(x, collapse=" ")
out
}))a_d <- ggplot(d1_lf, aes(x=scores, y=item_f)) +
geom_vline(xintercept=seq(0, 3, 1), size=1.5, color="gray57") +
geom_jitter(alpha=.3, size=2.5, width=.1, height=0.28, pch=21, color="darkred", fill="white") +
#geom_hline(yintercept=seq(1, 15, 1)) +
geom_errorbar(data=d1_summs, aes(x=m, y=item,
xmin=m-ci95, xmax=m+ci95), color="darkred",
inherit.aes=FALSE, size=2.2, width=.4) +
geom_point(data=d1_summs, aes(x=m, y=item), pch=21,
fill="white", color="black", size=6, alpha=.9) +
ga + ap + theme(axis.text.y =element_text(size=11)) +
theme(axis.text.x =element_text(size=15)) +
theme(axis.title.x =element_text(size=18)) +
ylab("") +
labs(title="Depression symptoms \nat baseline") + tp +
xlab("ratings") +
scale_y_discrete(labels=b_labs)#d1_summs # Spot check that life meaningless and life worthless are low means , easily upset in middle, couldnt get going at top a_d
Change
change_for_p <- data.frame("ID"=d2_short$ID, change_df[grep("change", names(change_df))])ch_lf <- change_for_p %>%
pivot_longer(cols = change_no.positive.feelings:change_no.initiative,
names_to = "item",
values_to = "scores"
)
ch_lf <- ch_lf[-c(which(ch_lf$scores %% 1 != 0)), ]
#which(ch_lf$scores %% 1 !=0)
ch_summs <- ch_lf %>% group_by(item) %>% summarize(m=mean(scores))
ch_sds <- ch_lf %>% group_by(item) %>%
summarize(ci95=qnorm(.975)*sd(scores)/sqrt(length(unique(bx_df$subj_idx))))
if (all(ch_summs$item==ch_sds$item)) ch_summs$ci95 <- ch_sds$ci95
ch_summs_arr <-
ch_summs %>% arrange(m)
ch_summs_arr$item_f <- factor(ch_summs_arr$item, levels=rev(ch_summs_arr$item))ch_lf$item_f <- factor(ch_lf$item, levels=rev(c(ch_summs_arr$item)))relab <- unlist(map(strsplit(levels(ch_lf$item_f), "change_"), 2))
sx_change_labs <- unlist(lapply(strsplit(unique(relab), "[.]"), function(x) {
out <- paste0(x, collapse=" ")
out
}))
sx_change_labs## [1] "life meaningless" "life not worthwhile"
## [3] "no hope" "down hearted and blue"
## [5] "no positive feelings" "lost interest"
## [7] "worthless" "not worth much as person"
## [9] "easily upset" "no initiative"
## [11] "no enthusiasm" "no enjoyment"
## [13] "sad and depressed" "couldnt get going"b_d <- ggplot(ch_lf, aes(x=scores, y=item_f)) +
geom_vline(xintercept=seq(-3, 3, 1), size=1.5, color="gray57") +
geom_jitter(alpha=.3, size=2.5, width=.1, height=0.28, pch=21, color="red", fill="white") +
#geom_hline(yintercept=seq(1, 15, 1)) +
geom_errorbar(data=ch_summs_arr, aes(x=m, y=item,
xmin=m-ci95, xmax=m+ci95), color="red",
inherit.aes=FALSE, size=2.2, width=.4) +
geom_point(data=ch_summs_arr, aes(x=m, y=item), pch=21,
fill="white", color="black", size=6, alpha=.8) +
ga + ap + theme(axis.text.y =element_text(size=11)) +
theme(axis.title.x =element_text(size=18)) +
theme(axis.text.x =element_text(size=15)) +
ylab("") +
labs(title="Change in depression symptoms \n(post - pre)") + tp +
xlab(TeX("$\\Delta$ ratings")) +
#scale_x_discrete(labels=bn) +
scale_y_discrete(labels=sx_change_labs)# Spot check
# names(d2_short)==names(d1_short)
# data.frame(colMeans(d2_short-d1_short)[1:14])b_d
ch_lf## # A tibble: 1,296 × 4
## ID item scores item_f
## <int> <chr> <dbl> <fct>
## 1 1 change_no.positive.feelings -1 change_no.positive.feelings
## 2 1 change_couldnt.get.going -1 change_couldnt.get.going
## 3 1 change_easily.upset -2 change_easily.upset
## 4 1 change_sad.and.depressed -1 change_sad.and.depressed
## 5 1 change_lost.interest -1 change_lost.interest
## 6 1 change_not.worth.much.as.person 1 change_not.worth.much.as.person
## 7 1 change_life.not.worthwhile 0 change_life.not.worthwhile
## 8 1 change_no.enjoyment 0 change_no.enjoyment
## 9 1 change_down.hearted.and.blue 0 change_down.hearted.and.blue
## 10 1 change_no.enthusiasm -1 change_no.enthusiasm
## # … with 1,286 more rowsunique_items <- unique(d1_lf$item)
d1m_s <- d1m %>% filter(ID %in% d1_lf$ID)
all_res <- foreach (i = 1:length(unique_items)) %do% {
this_item <- d1_lf %>% filter(item==unique_items[i])
d1m_s_tmp <- d1m_s %>% filter(ID %in% unique(this_item$ID))
#print(this_item)
if (all(this_item$ID==d1m_s_tmp$ID)) {
res <- cor.test(this_item$scores, d1m_s_tmp$v_val_ctr_)
out <- data.table("item"=unique(this_item$item),
"r"=as.numeric(res$estimate),
"ci95_low"=as.numeric(res$conf.int[1]),
"ci95_high"=as.numeric(res$conf.int[2]),
"p"=as.numeric(res$p.value))
} else {
}
out
} %>% bind_rows()bonferroni <- .05/length(unique(d1_lf$item))
all_res$bonferroni <- (all_res$p < bonferroni)*1
all_res[which(all_res$p < .05), "signif"] <- "*"
all_res[which(all_res$p < .01), "signif"] <- "**"
all_res[which(all_res$p < .005), "signif"] <- "***"
all_res[which(all_res$p < .001), "signif"] <- "****"
all_res[which(all_res$p < .0005), "signif"] <- "*****"
all_res[which(is.na(all_res$signif)), "signif"] <- " "all_res_arr <- all_res %>% arrange(r)
all_res_arr$item <- factor(all_res_arr$item, levels=rev(c(all_res_arr$item)))#all_res_arr %>% select(item, bonferroni, signif)#ch_lf$item_f <- factor(ch_lf$item, levels=rev(c(ch_summs_arr$item)))
#relab <- unlist(map(strsplit(levels(ch_lf$item_f), "change_"), 2))
r_b_labs <- as.character(levels(all_res_arr$item))
rb_labs <- unlist(lapply(strsplit(r_b_labs, "[.]"), function(x) {
out <- paste0(x, collapse=" ")
out
}))a_r <- ggplot(all_res_arr, aes(x=r, y=item)) +
geom_vline(xintercept = seq(-1, .2, .1), color="gray57", linetype="dotted") +
geom_vline(xintercept = c(-1, 0), color="black", size=1.4) +
geom_errorbar(data=all_res_arr, size=1.4, width=.19, color="gray42",
aes(x=r, y=item, xmin=ci95_low, xmax=ci95_high)) +
#geom_hline(yintercept=seq(1, 15, 1)) +
geom_point(size=6, fill="gray60", pch=21, alpha=.9) +
xlim(-1.15, .22) +
geom_text(x=-1.1, y=14, size=6, label="X*****") +
geom_text(x=-1.1, y=13, size=6, label="X*****") +
geom_text(x=-1.1, y=12, size=6, label="X ***") +
geom_text(x=-1.1, y=11, size=6, label="X ***") +
geom_text(x=-1.1, y=10, size=6, label="X ***") +
geom_text(x=-1.1, y=9, size=6, label="X ***") +
geom_text(x=-1.1, y=8, size=6, label=" **") +
geom_text(x=-1.1, y=7, size=6, label=" **") +
geom_text(x=-1.1, y=6, size=6, label=" *") +
geom_text(x=-1.1, y=5, size=6, label=" *") +
geom_text(x=-1.1, y=4, size=6, label=" *") +
geom_text(x=-1.1, y=3, size=6, label=" *") +
ga + theme(axis.text.y =element_text(size=11), axis.text.x=element_text(size=15),
axis.title = element_text(size=20)) +
ylab("") +
labs(title="Correlations of baseline depression symptoms with \n positive > negative drift rate estimates") +
theme(plot.title = element_text(size = 15, face='bold', hjust = .5)) +
xlab("Pearson's r (with 95% CIs)") +
scale_y_discrete(labels=rb_labs)
a_r
unique_ch_items <- unique(ch_lf$item)
d1m_s_ch <- d1m %>% filter(ID %in% ch_lf$ID)
#d1m_s_ch
all_res_ch <- foreach (i = 1:length(unique_ch_items)) %do% {
this_ch_item <- ch_lf %>% filter(item==unique_ch_items[i])
d1m_s_ch_tmp <- d1m_s_ch %>% filter(ID %in% unique(this_ch_item$ID))
#print(this_item)
if (all(this_ch_item$ID==d1m_s_ch_tmp$ID)) {
res <- cor.test(this_ch_item$scores, d1m_s_ch_tmp$v_val_ctr.sess_ctr_)
out <- data.table("item"=unique(this_ch_item$item),
"r"=as.numeric(res$estimate),
"ci95_low"=as.numeric(res$conf.int[1]),
"ci95_high"=as.numeric(res$conf.int[2]),
"p"=as.numeric(res$p.value))
} else {
}
out
} %>% bind_rows()bonferroni_ch <- .05/length(unique(ch_lf$item))
all_res_ch$bonferroni_ch <- (all_res_ch$p < bonferroni_ch)*1
all_res_ch[which(all_res_ch$p < .05), "signif"] <- "*"
all_res_ch[which(all_res_ch$p < .01), "signif"] <- "**"
all_res_ch[which(all_res_ch$p < .005), "signif"] <- "***"
all_res_ch[which(all_res_ch$p < .001), "signif"] <- "****"
all_res_ch[which(is.na(all_res_ch$signif)), "signif"] <- " "
#all_res_ch
arc_arr <- all_res_ch %>% arrange(r)
#arc_arr
arc_arr$item_f <- factor(arc_arr$item, levels=rev(as.character(arc_arr$item)))
#arc_arr$item_f <- factor(arc_arr$item, levels=as.character(arc_arr$item))
#arc_arr$item_fr_ch_labs <- as.character(levels(arc_arr$item_f))
rch_labs <- unlist(lapply(strsplit(r_ch_labs, "[.]"), function(x) {
out <- paste0(x, collapse=" ")
out <- unlist(map(strsplit(out, "change_"), 2))
out
}))# all_res_ch # Spot checked b_r <- ggplot(arc_arr, aes(x=r, y=item_f)) +
geom_text(x=-1.1, y=14, size=6, label="X ***") +
geom_text(x=-1.1, y=13, size=6, label="X ***") +
geom_text(x=-1.1, y=12, size=6, label=" ***") +
geom_text(x=-1.1, y=11, size=6, label=" **") +
# geom_text(x=-.35, y=10, size=6, label=" **") +
# geom_text(x=-.35, y=9, size=6, label=" *") +
geom_vline(xintercept = seq(-1, .2, .1), color="gray57", linetype="dotted") +
geom_vline(xintercept = c(-1, 0), color="black", size=1.4) +
xlim(-1.15, .22) +
geom_errorbar(data=arc_arr, size=1.4, width=.19, color="gray75",
aes(x=r, y=item_f, xmin=ci95_low, xmax=ci95_high)) +
geom_point(size=6, fill="gray85", pch=21, alpha=.9) +
ga +
theme(axis.text.y =element_text(size=11), axis.text.x=element_text(size=15),
axis.title = element_text(size=20)) +
ylab("") +
labs(title="Correlations of change in depression symptoms with \n positive > negative drift rate * time estimates") +
theme(plot.title = element_text(size = 15, face='bold', hjust = .5)) +
xlab("Pearson's r (with 95% CIs)") +
scale_y_discrete(labels=rch_labs)
b_r
full_sxs_p <- a_d / b_d/
a_r / b_r
#full_sxs_p#ggsave("../paper/figs/full_sxs.png", full_sxs_p, dpi=700, width=8, height=16)Get Cronbach’s alpha
psych::alpha(dd1[2:15])##
## Reliability analysis
## Call: psych::alpha(x = dd1[2:15])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.93 0.93 0.95 0.5 14 0.011 0.68 0.55 0.47
##
## lower alpha upper 95% confidence boundaries
## 0.91 0.93 0.95
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## DASS1.3 0.93 0.93 0.95 0.50 13 0.011 0.012 0.48
## DASS1.5 0.93 0.93 0.95 0.51 13 0.011 0.012 0.49
## DASS1.10 0.92 0.92 0.95 0.49 12 0.012 0.012 0.46
## DASS1.13 0.93 0.93 0.95 0.51 13 0.011 0.011 0.49
## DASS1.16 0.92 0.93 0.95 0.49 13 0.012 0.012 0.47
## DASS1.17 0.93 0.93 0.94 0.50 13 0.011 0.011 0.48
## DASS1.21 0.93 0.93 0.94 0.50 13 0.011 0.010 0.47
## DASS1.24 0.92 0.92 0.94 0.49 12 0.012 0.012 0.47
## DASS1.26 0.92 0.93 0.94 0.49 13 0.012 0.012 0.47
## DASS1.31 0.92 0.93 0.95 0.50 13 0.011 0.011 0.47
## DASS1.34 0.92 0.93 0.94 0.49 13 0.011 0.011 0.47
## DASS1.37 0.92 0.93 0.95 0.49 12 0.012 0.011 0.47
## DASS1.38 0.92 0.93 0.94 0.49 12 0.012 0.011 0.47
## DASS1.42 0.93 0.93 0.95 0.51 13 0.011 0.011 0.49
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## DASS1.3 95 0.70 0.70 0.67 0.64 0.66 0.78
## DASS1.5 95 0.66 0.64 0.61 0.59 1.20 0.86
## DASS1.10 95 0.80 0.80 0.79 0.76 0.63 0.85
## DASS1.13 95 0.65 0.65 0.62 0.58 1.03 0.82
## DASS1.16 95 0.75 0.75 0.72 0.70 0.48 0.67
## DASS1.17 95 0.69 0.70 0.69 0.64 0.57 0.75
## DASS1.21 95 0.70 0.72 0.71 0.66 0.32 0.61
## DASS1.24 95 0.80 0.80 0.79 0.76 0.78 0.79
## DASS1.26 95 0.76 0.75 0.74 0.71 0.88 0.78
## DASS1.31 95 0.73 0.73 0.71 0.68 0.72 0.79
## DASS1.34 95 0.74 0.75 0.74 0.69 0.49 0.74
## DASS1.37 95 0.78 0.79 0.77 0.74 0.43 0.69
## DASS1.38 95 0.75 0.77 0.77 0.72 0.25 0.56
## DASS1.42 95 0.67 0.65 0.62 0.60 1.13 0.88
##
## Non missing response frequency for each item
## 0 0.384615384615385 1 1.61538461538462 1.92307692307692 2
## DASS1.3 0.48 0.00 0.41 0.00 0.00 0.06
## DASS1.5 0.19 0.00 0.52 0.00 0.00 0.20
## DASS1.10 0.57 0.00 0.27 0.00 0.00 0.12
## DASS1.13 0.25 0.01 0.52 0.00 0.00 0.16
## DASS1.16 0.60 0.00 0.33 0.00 0.00 0.06
## DASS1.17 0.57 0.00 0.32 0.00 0.00 0.09
## DASS1.21 0.75 0.00 0.20 0.00 0.00 0.04
## DASS1.24 0.42 0.00 0.40 0.00 0.00 0.16
## DASS1.26 0.33 0.00 0.51 0.00 0.00 0.13
## DASS1.31 0.46 0.00 0.39 0.00 0.00 0.12
## DASS1.34 0.62 0.00 0.29 0.00 0.01 0.04
## DASS1.37 0.66 0.00 0.26 0.00 0.00 0.05
## DASS1.38 0.80 0.00 0.16 0.00 0.00 0.03
## DASS1.42 0.24 0.00 0.46 0.01 0.00 0.20
## 3 miss
## DASS1.3 0.04 0
## DASS1.5 0.09 0
## DASS1.10 0.04 0
## DASS1.13 0.06 0
## DASS1.16 0.01 0
## DASS1.17 0.02 0
## DASS1.21 0.01 0
## DASS1.24 0.02 0
## DASS1.26 0.04 0
## DASS1.31 0.03 0
## DASS1.34 0.03 0
## DASS1.37 0.02 0
## DASS1.38 0.01 0
## DASS1.42 0.08 0psych::alpha(dd2[2:15])##
## Reliability analysis
## Call: psych::alpha(x = dd2[2:15])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.93 0.94 0.96 0.51 15 0.01 0.31 0.38 0.52
##
## lower alpha upper 95% confidence boundaries
## 0.91 0.93 0.95
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## DASS2.3 0.92 0.93 0.96 0.51 14 0.011 0.022 0.51
## DASS2.5 0.92 0.93 0.96 0.52 14 0.011 0.021 0.52
## DASS2.10 0.92 0.93 0.96 0.51 13 0.011 0.020 0.52
## DASS2.13 0.92 0.93 0.95 0.50 13 0.012 0.021 0.51
## DASS2.16 0.92 0.93 0.96 0.50 13 0.011 0.019 0.51
## DASS2.17 0.93 0.94 0.96 0.53 15 0.010 0.018 0.53
## DASS2.21 0.92 0.93 0.96 0.50 13 0.011 0.021 0.50
## DASS2.24 0.92 0.93 0.96 0.51 13 0.011 0.019 0.51
## DASS2.26 0.92 0.93 0.96 0.51 13 0.011 0.022 0.52
## DASS2.31 0.92 0.93 0.96 0.51 14 0.011 0.020 0.51
## DASS2.34 0.93 0.94 0.96 0.54 15 0.010 0.017 0.53
## DASS2.37 0.92 0.93 0.95 0.50 13 0.011 0.020 0.51
## DASS2.38 0.92 0.93 0.96 0.51 14 0.011 0.022 0.52
## DASS2.42 0.92 0.93 0.96 0.52 14 0.011 0.021 0.53
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## DASS2.3 94 0.76 0.74 0.72 0.71 0.35 0.56
## DASS2.5 94 0.75 0.70 0.68 0.67 0.70 0.76
## DASS2.10 94 0.74 0.76 0.75 0.70 0.20 0.50
## DASS2.13 94 0.82 0.80 0.79 0.78 0.50 0.63
## DASS2.16 94 0.81 0.84 0.84 0.78 0.18 0.41
## DASS2.17 94 0.54 0.58 0.55 0.49 0.16 0.40
## DASS2.21 94 0.79 0.82 0.81 0.77 0.14 0.35
## DASS2.24 94 0.76 0.77 0.76 0.71 0.29 0.58
## DASS2.26 94 0.79 0.76 0.75 0.74 0.59 0.68
## DASS2.31 94 0.75 0.76 0.74 0.70 0.22 0.47
## DASS2.34 94 0.52 0.55 0.51 0.46 0.13 0.40
## DASS2.37 94 0.79 0.81 0.81 0.76 0.16 0.40
## DASS2.38 94 0.70 0.74 0.72 0.66 0.11 0.34
## DASS2.42 94 0.77 0.72 0.70 0.70 0.67 0.77
##
## Non missing response frequency for each item
## 0 0.636363636363636 1 2 3 miss
## DASS2.3 0.69 0.00 0.27 0.04 0.00 0
## DASS2.5 0.44 0.01 0.46 0.05 0.04 0
## DASS2.10 0.83 0.00 0.15 0.01 0.01 0
## DASS2.13 0.57 0.01 0.34 0.07 0.00 0
## DASS2.16 0.83 0.00 0.16 0.01 0.00 0
## DASS2.17 0.85 0.00 0.14 0.01 0.00 0
## DASS2.21 0.86 0.00 0.14 0.00 0.00 0
## DASS2.24 0.77 0.00 0.19 0.03 0.01 0
## DASS2.26 0.52 0.00 0.37 0.11 0.00 0
## DASS2.31 0.80 0.00 0.18 0.02 0.00 0
## DASS2.34 0.88 0.01 0.09 0.02 0.00 0
## DASS2.37 0.85 0.00 0.14 0.01 0.00 0
## DASS2.38 0.90 0.00 0.09 0.01 0.00 0
## DASS2.42 0.48 0.00 0.40 0.09 0.03 0Continuing to prep for regressions…
FFMQ
iq_c_m <- iq_c[grep("FFMQ", names(iq_c))]
# Reverse code items
rev_items <- c(3, 5, 8, 10, 12, 13, 14, 16, 17, 18, 22, 23, 25, 28, 30, 34, 35, 38, 39)
# Check
#sort(c(3, 10, 14, 17, 25, 30, 35, 39, 5, 8, 13, 18, 23, 28, 34, 38, 12, 16, 22))==rev_items
aware <- c(5, 8, 13, 18, 23, 28, 34, 38)
non_judge <- c(3, 10, 14, 17, 25, 30, 35, 39)
non_react <- c(4, 9, 19, 21, 24, 29, 33)
observe <- c(1, 6, 11, 15, 20, 26, 31, 36)Time 1
m1 <- iq_c_m[grep("FFMQ1.", names(iq_c_m))]
# 1.31% of baseline data missing
sum(length(which(m1==999)), length(which(is.na(m1))))/(dim(m1)[1]*dim(m1)[2])*100## [1] 1.308761m1[m1==999] <- NA
m1[paste0("FFMQ1.", rev_items)] <- 6-m1[paste0("FFMQ1.", rev_items)]
m1 <- data.frame("ID"=iq_c$ID.ffmq, m1)
#m1[which(is.na(rowSums(m1))), ]
# Delete 50 who's also missing ffmq data
m1 <- m1[!m1$ID==50, ]
length(which(is.na(m1)))/(dim(m1)[1]*dim(m1)[2])*100## [1] 0.2631579# For robustness check, casewise delete pts with any data missing
m1_casewise_delete <- m1[-which(is.na(rowSums(m1))), ]
# Replace any NAs with the mean of the remainding items
for (r in 1:nrow(m1)) {
if (any(is.na(m1[r, 2:40]))) {
this_row <- as.numeric(unlist(m1[r, 2:40])) # Exclude the ID
# Mean to replace NAs with
mean_non_na <- mean(this_row[which(!is.na(this_row))])
this_row[is.na(as.numeric(unlist(this_row)))] <- mean_non_na
# Put in df
m1[r, 2:40] <- this_row
}
}
#m1[which(is.na(rowSums(m1))), ]
#which(is.na(m1))
m1$aware_sum <- rowSums(m1[2:40][aware])
m1$nonjudge_sum <- rowSums(m1[2:40][non_judge])
m1$nonreact_sum <- rowSums(m1[2:40][non_react])
m1$observe_sum <- rowSums(m1[2:40][observe])
m1_casewise_delete$aware_sum <- rowSums(m1_casewise_delete[2:40][aware])
m1_casewise_delete$nonjudge_sum <- rowSums(m1_casewise_delete[2:40][non_judge])
m1_casewise_delete$nonreact_sum <- rowSums(m1_casewise_delete[2:40][non_react])
m1_casewise_delete$observe_sum <- rowSums(m1_casewise_delete[2:40][observe])
cor.test(m1$observe_sum, m1$aware_sum)##
## Pearson's product-moment correlation
##
## data: m1$observe_sum and m1$aware_sum
## t = 0.60707, df = 93, p-value = 0.5453
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1404957 0.2610637
## sample estimates:
## cor
## 0.06282629cor.test(m1$observe_sum, m1$nonjudge_sum)##
## Pearson's product-moment correlation
##
## data: m1$observe_sum and m1$nonjudge_sum
## t = 0.35017, df = 93, p-value = 0.727
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1664736 0.2361032
## sample estimates:
## cor
## 0.03628692cor.test(m1$observe_sum, m1$nonreact_sum)##
## Pearson's product-moment correlation
##
## data: m1$observe_sum and m1$nonreact_sum
## t = 2.4114, df = 93, p-value = 0.01786
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.04314921 0.42342407
## sample estimates:
## cor
## 0.2425826cor.test(m1$aware_sum, m1$nonreact_sum)##
## Pearson's product-moment correlation
##
## data: m1$aware_sum and m1$nonreact_sum
## t = 3.6369, df = 93, p-value = 0.0004525
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1629167 0.5176046
## sample estimates:
## cor
## 0.3528732mean(m1$aware_sum); sd(m1$aware_sum)## [1] 24.04737## [1] 5.984788mean(m1$nonjudge_sum); sd(m1$nonjudge_sum)## [1] 26.30471## [1] 7.658276mean(m1$nonreact_sum); sd(m1$nonreact_sum)## [1] 18.83914## [1] 3.837917psych::alpha(m1[2:40][aware])##
## Reliability analysis
## Call: psych::alpha(x = m1[2:40][aware])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.88 0.88 0.9 0.48 7.5 0.018 3 0.75 0.52
##
## lower alpha upper 95% confidence boundaries
## 0.85 0.88 0.92
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## FFMQ1.5 0.86 0.86 0.87 0.47 6.3 0.022 0.020 0.52
## FFMQ1.8 0.86 0.86 0.88 0.46 6.0 0.023 0.028 0.52
## FFMQ1.13 0.86 0.86 0.87 0.47 6.3 0.022 0.021 0.52
## FFMQ1.18 0.88 0.88 0.90 0.51 7.3 0.019 0.027 0.54
## FFMQ1.23 0.87 0.87 0.89 0.50 6.9 0.020 0.030 0.52
## FFMQ1.28 0.87 0.87 0.89 0.48 6.5 0.021 0.031 0.53
## FFMQ1.34 0.88 0.88 0.90 0.52 7.6 0.018 0.020 0.53
## FFMQ1.38 0.86 0.85 0.87 0.46 5.9 0.023 0.029 0.46
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## FFMQ1.5 95 0.79 0.78 0.78 0.71 2.8 1.08
## FFMQ1.8 95 0.82 0.82 0.79 0.76 3.3 1.01
## FFMQ1.13 95 0.80 0.78 0.77 0.71 2.9 1.15
## FFMQ1.18 95 0.65 0.65 0.57 0.54 3.0 0.94
## FFMQ1.23 95 0.70 0.70 0.64 0.59 3.0 1.04
## FFMQ1.28 95 0.74 0.75 0.70 0.65 3.1 0.96
## FFMQ1.34 95 0.59 0.61 0.54 0.48 3.1 0.91
## FFMQ1.38 95 0.84 0.84 0.83 0.78 3.0 0.95
##
## Non missing response frequency for each item
## 1 2 2.71052631578947 2.78947368421053 3 4 5 miss
## FFMQ1.5 0.14 0.26 0.00 0.00 0.33 0.23 0.04 0
## FFMQ1.8 0.01 0.24 0.00 0.00 0.32 0.31 0.13 0
## FFMQ1.13 0.16 0.19 0.00 0.00 0.34 0.25 0.06 0
## FFMQ1.18 0.04 0.27 0.01 0.00 0.40 0.22 0.05 0
## FFMQ1.23 0.05 0.29 0.00 0.01 0.34 0.22 0.08 0
## FFMQ1.28 0.02 0.28 0.00 0.00 0.36 0.26 0.07 0
## FFMQ1.34 0.02 0.23 0.00 0.00 0.45 0.22 0.07 0
## FFMQ1.38 0.05 0.25 0.00 0.00 0.42 0.22 0.05 0psych::alpha(m1[2:40][non_judge])##
## Reliability analysis
## Call: psych::alpha(x = m1[2:40][non_judge])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.94 0.94 0.94 0.67 16 0.0091 3.3 0.96 0.67
##
## lower alpha upper 95% confidence boundaries
## 0.92 0.94 0.96
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## FFMQ1.3 0.94 0.94 0.94 0.69 16 0.0096 0.0035 0.68
## FFMQ1.10 0.93 0.94 0.94 0.67 14 0.0103 0.0052 0.67
## FFMQ1.14 0.94 0.94 0.93 0.68 15 0.0101 0.0040 0.67
## FFMQ1.17 0.93 0.93 0.93 0.66 13 0.0110 0.0035 0.65
## FFMQ1.25 0.93 0.93 0.93 0.66 13 0.0110 0.0041 0.65
## FFMQ1.30 0.93 0.93 0.93 0.67 14 0.0106 0.0042 0.66
## FFMQ1.35 0.93 0.93 0.93 0.67 14 0.0104 0.0036 0.66
## FFMQ1.39 0.94 0.94 0.94 0.68 15 0.0101 0.0053 0.67
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## FFMQ1.3 95 0.78 0.78 0.74 0.72 3.1 1.1
## FFMQ1.10 95 0.84 0.84 0.81 0.78 3.1 1.1
## FFMQ1.14 95 0.83 0.83 0.80 0.77 3.6 1.1
## FFMQ1.17 95 0.89 0.89 0.88 0.85 3.1 1.2
## FFMQ1.25 95 0.88 0.89 0.87 0.85 3.1 1.1
## FFMQ1.30 95 0.86 0.86 0.84 0.81 3.5 1.1
## FFMQ1.35 95 0.85 0.84 0.83 0.79 3.4 1.2
## FFMQ1.39 95 0.83 0.83 0.79 0.77 3.4 1.2
##
## Non missing response frequency for each item
## 1 2 2.94736842105263 3 4 5 miss
## FFMQ1.3 0.06 0.25 0.00 0.34 0.23 0.12 0
## FFMQ1.10 0.07 0.24 0.00 0.36 0.21 0.12 0
## FFMQ1.14 0.02 0.15 0.00 0.36 0.16 0.32 0
## FFMQ1.17 0.04 0.35 0.00 0.28 0.17 0.16 0
## FFMQ1.25 0.03 0.29 0.00 0.35 0.18 0.15 0
## FFMQ1.30 0.02 0.16 0.01 0.32 0.25 0.24 0
## FFMQ1.35 0.04 0.19 0.00 0.31 0.22 0.24 0
## FFMQ1.39 0.03 0.25 0.00 0.24 0.20 0.27 0psych::alpha(m1[2:40][non_react])##
## Reliability analysis
## Call: psych::alpha(x = m1[2:40][non_react])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.77 0.78 0.78 0.33 3.5 0.036 2.7 0.55 0.32
##
## lower alpha upper 95% confidence boundaries
## 0.7 0.77 0.84
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## FFMQ1.4 0.77 0.77 0.77 0.36 3.4 0.037 0.018 0.39
## FFMQ1.9 0.72 0.73 0.73 0.31 2.7 0.045 0.023 0.28
## FFMQ1.19 0.74 0.75 0.73 0.33 2.9 0.041 0.015 0.30
## FFMQ1.21 0.75 0.76 0.76 0.34 3.1 0.040 0.021 0.34
## FFMQ1.24 0.78 0.78 0.78 0.38 3.6 0.036 0.017 0.39
## FFMQ1.29 0.71 0.71 0.71 0.29 2.5 0.047 0.015 0.28
## FFMQ1.33 0.73 0.73 0.73 0.31 2.7 0.043 0.018 0.32
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## FFMQ1.4 95 0.58 0.55 0.44 0.38 2.8 0.96
## FFMQ1.9 95 0.73 0.73 0.67 0.60 2.7 0.87
## FFMQ1.19 95 0.65 0.66 0.60 0.50 2.7 0.81
## FFMQ1.21 95 0.62 0.62 0.53 0.46 3.0 0.85
## FFMQ1.24 95 0.50 0.50 0.36 0.31 2.7 0.84
## FFMQ1.29 95 0.78 0.78 0.77 0.66 2.6 0.85
## FFMQ1.33 95 0.70 0.72 0.68 0.58 2.4 0.72
##
## Non missing response frequency for each item
## 1 2 3 3.2972972972973 3.68421052631579 3.73684210526316 4
## FFMQ1.4 0.08 0.33 0.37 0.00 0.00 0.00 0.19
## FFMQ1.9 0.03 0.47 0.31 0.00 0.00 0.00 0.17
## FFMQ1.19 0.05 0.37 0.42 0.01 0.01 0.00 0.13
## FFMQ1.21 0.00 0.28 0.44 0.00 0.00 0.00 0.22
## FFMQ1.24 0.06 0.37 0.39 0.00 0.00 0.01 0.17
## FFMQ1.29 0.07 0.38 0.40 0.00 0.00 0.00 0.14
## FFMQ1.33 0.07 0.55 0.32 0.00 0.00 0.00 0.06
## 5 miss
## FFMQ1.4 0.03 0
## FFMQ1.9 0.02 0
## FFMQ1.19 0.01 0
## FFMQ1.21 0.05 0
## FFMQ1.24 0.00 0
## FFMQ1.29 0.01 0
## FFMQ1.33 0.00 0psych::alpha(m1[2:40][observe])##
## Reliability analysis
## Call: psych::alpha(x = m1[2:40][observe])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.82 0.82 0.84 0.37 4.7 0.028 3.2 0.69 0.35
##
## lower alpha upper 95% confidence boundaries
## 0.76 0.82 0.88
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## FFMQ1.1 0.80 0.80 0.82 0.37 4.1 0.032 0.019 0.35
## FFMQ1.6 0.78 0.79 0.80 0.35 3.7 0.034 0.016 0.33
## FFMQ1.11 0.82 0.82 0.82 0.40 4.6 0.028 0.013 0.41
## FFMQ1.15 0.79 0.79 0.80 0.35 3.8 0.034 0.016 0.35
## FFMQ1.20 0.80 0.81 0.81 0.38 4.2 0.031 0.013 0.35
## FFMQ1.26 0.80 0.80 0.81 0.37 4.1 0.032 0.017 0.35
## FFMQ1.31 0.79 0.80 0.81 0.36 3.9 0.033 0.017 0.35
## FFMQ1.36 0.81 0.82 0.82 0.39 4.5 0.030 0.016 0.41
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## FFMQ1.1 95 0.67 0.67 0.60 0.54 2.5 1.03
## FFMQ1.6 95 0.75 0.75 0.72 0.64 2.8 1.02
## FFMQ1.11 95 0.59 0.56 0.48 0.41 3.0 1.22
## FFMQ1.15 95 0.74 0.75 0.71 0.64 3.2 0.92
## FFMQ1.20 95 0.63 0.64 0.59 0.50 3.4 0.98
## FFMQ1.26 95 0.67 0.68 0.62 0.55 3.8 0.94
## FFMQ1.31 95 0.71 0.71 0.67 0.59 3.8 1.07
## FFMQ1.36 95 0.61 0.59 0.52 0.47 3.3 1.04
##
## Non missing response frequency for each item
## 1 2 3 3.2972972972973 4 5 miss
## FFMQ1.1 0.19 0.34 0.28 0.00 0.18 0.01 0
## FFMQ1.6 0.13 0.24 0.37 0.00 0.24 0.02 0
## FFMQ1.11 0.14 0.23 0.23 0.00 0.31 0.09 0
## FFMQ1.15 0.05 0.13 0.47 0.00 0.28 0.06 0
## FFMQ1.20 0.04 0.11 0.42 0.00 0.31 0.13 0
## FFMQ1.26 0.02 0.06 0.24 0.01 0.43 0.23 0
## FFMQ1.31 0.03 0.07 0.25 0.00 0.32 0.33 0
## FFMQ1.36 0.05 0.17 0.36 0.00 0.31 0.12 0# Time 2
# 2.24 % of final data missing
m2 <- iq_c_m[grep("FFMQ2.", names(iq_c_m))]
sum(length(which(m2==999)), length(which(is.na(m2))))/(dim(m2)[1]*dim(m2)[2])*100## [1] 2.24359m2[m2==999] <- NA
m2[paste0("FFMQ2.", rev_items)] <- 6-m2[paste0("FFMQ2.", rev_items)]
m2 <- data.frame("ID"=iq_c$ID.ffmq, m2)
m2[which(is.na(rowSums(m2))), ]## ID FFMQ2.1 FFMQ2.2 FFMQ2.3 FFMQ2.4 FFMQ2.5 FFMQ2.6 FFMQ2.7 FFMQ2.8 FFMQ2.9
## 30 34 3 3 5 3 3 1 2 3 3
## 38 44 4 3 3 3 4 NA 3 4 3
## 48 56 NA NA NA NA NA NA NA NA NA
## 50 58 NA NA NA NA NA NA NA NA NA
## 86 101 4 5 3 3 4 4 5 3 3
## 92 109 5 5 5 1 5 5 5 5 4
## 96 113 4 4 5 NA 2 4 3 3 4
## FFMQ2.10 FFMQ2.11 FFMQ2.12 FFMQ2.13 FFMQ2.14 FFMQ2.15 FFMQ2.16 FFMQ2.17
## 30 5 5 4 3 5 2 3 5
## 38 3 4 3 4 3 4 NA 3
## 48 NA NA NA NA NA NA NA NA
## 50 NA NA NA NA NA NA NA NA
## 86 4 4 5 3 4 5 5 4
## 92 5 4 5 5 5 5 5 5
## 96 5 5 4 3 5 5 3 4
## FFMQ2.18 FFMQ2.19 FFMQ2.20 FFMQ2.21 FFMQ2.22 FFMQ2.23 FFMQ2.24 FFMQ2.25
## 30 3 4 2 3 2 2 3 5
## 38 4 3 4 3 4 4 3 3
## 48 NA NA NA NA NA NA NA NA
## 50 NA NA NA NA NA NA NA NA
## 86 3 3 5 4 5 3 3 3
## 92 5 4 5 NA 5 5 4 5
## 96 3 4 5 4 4 4 4 5
## FFMQ2.26 FFMQ2.27 FFMQ2.28 FFMQ2.29 FFMQ2.30 FFMQ2.31 FFMQ2.32 FFMQ2.33
## 30 2 2 3 3 5 NA 1 3
## 38 3 3 4 3 3 4 3 3
## 48 NA NA NA NA NA NA NA NA
## 50 NA NA NA NA NA NA NA NA
## 86 5 5 3 3 5 5 5 2
## 92 5 5 5 4 5 5 5 4
## 96 5 3 3 4 5 5 4 5
## FFMQ2.34 FFMQ2.35 FFMQ2.36 FFMQ2.37 FFMQ2.38 FFMQ2.39
## 30 3 5 4 2 2 5
## 38 4 4 3 3 4 3
## 48 NA NA NA NA NA NA
## 50 NA NA NA NA NA NA
## 86 3 NA 5 5 3 4
## 92 5 5 5 5 5 5
## 96 3 5 5 3 3 5ids_to_del <- c(56, 58)
# Delete 50 who's also missing ffmq data
m2 <- m2 %>% filter(!ID %in% ids_to_del)
length(which(is.na(m2)))/(dim(m2)[1]*dim(m2)[2])*100## [1] 0.1595745# For robustness check, casewise delete pts with any data missing
m2_casewise_delete <- m2[-which(is.na(rowSums(m2))), ]
# Replace any NAs with the mean of the remainding items
for (r in 1:nrow(m2)) {
if (any(is.na(m2[r, 2:40]))) {
this_row <- as.numeric(unlist(m2[r, 2:40])) # Exclude the ID
# Mean to replace NAs with
mean_non_na <- mean(this_row[which(!is.na(this_row))])
this_row[is.na(as.numeric(unlist(this_row)))] <- mean_non_na
# Put in df
m2[r, 2:40] <- this_row
}
}
#which(is.na(m2))
m2$aware_sum <- rowSums(m2[2:40][aware])
m2$nonjudge_sum <- rowSums(m2[2:40][non_judge])
m2$nonreact_sum <- rowSums(m2[2:40][non_react])
m2$observe_sum <- rowSums(m2[2:40][observe])
m2_casewise_delete$aware_sum <- rowSums(m2_casewise_delete[2:40][aware])
m2_casewise_delete$nonjudge_sum <- rowSums(m2_casewise_delete[2:40][non_judge])
m2_casewise_delete$nonreact_sum <- rowSums(m2_casewise_delete[2:40][non_react])psych::alpha(m2[2:40][aware])##
## Reliability analysis
## Call: psych::alpha(x = m2[2:40][aware])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.9 0.9 0.91 0.53 9.1 0.016 3.5 0.63 0.54
##
## lower alpha upper 95% confidence boundaries
## 0.87 0.9 0.93
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## FFMQ2.5 0.89 0.89 0.90 0.54 8.2 0.018 0.013 0.56
## FFMQ2.8 0.88 0.88 0.90 0.52 7.6 0.019 0.016 0.55
## FFMQ2.13 0.89 0.89 0.89 0.53 7.9 0.018 0.012 0.54
## FFMQ2.18 0.89 0.89 0.91 0.54 8.3 0.018 0.016 0.56
## FFMQ2.23 0.89 0.89 0.90 0.53 8.0 0.017 0.014 0.55
## FFMQ2.28 0.89 0.89 0.90 0.55 8.5 0.017 0.012 0.55
## FFMQ2.34 0.88 0.88 0.90 0.52 7.6 0.019 0.014 0.54
## FFMQ2.38 0.88 0.88 0.89 0.51 7.2 0.019 0.015 0.50
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## FFMQ2.5 94 0.75 0.73 0.70 0.66 3.2 0.87
## FFMQ2.8 94 0.80 0.80 0.76 0.73 3.6 0.83
## FFMQ2.13 94 0.80 0.77 0.75 0.70 3.4 0.99
## FFMQ2.18 94 0.74 0.73 0.67 0.64 3.4 0.84
## FFMQ2.23 94 0.74 0.76 0.72 0.66 3.6 0.84
## FFMQ2.28 94 0.68 0.71 0.66 0.60 3.5 0.70
## FFMQ2.34 94 0.79 0.81 0.78 0.73 3.6 0.73
## FFMQ2.38 94 0.84 0.84 0.83 0.78 3.5 0.80
##
## Non missing response frequency for each item
## 1 2 3 4 5 miss
## FFMQ2.5 0.03 0.15 0.46 0.31 0.05 0
## FFMQ2.8 0.02 0.03 0.38 0.44 0.13 0
## FFMQ2.13 0.03 0.15 0.38 0.31 0.13 0
## FFMQ2.18 0.01 0.10 0.43 0.37 0.10 0
## FFMQ2.23 0.00 0.09 0.40 0.37 0.14 0
## FFMQ2.28 0.00 0.04 0.52 0.36 0.07 0
## FFMQ2.34 0.00 0.04 0.46 0.40 0.10 0
## FFMQ2.38 0.00 0.09 0.44 0.37 0.11 0psych::alpha(m2[2:40][non_judge])##
## Reliability analysis
## Call: psych::alpha(x = m2[2:40][non_judge])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.94 0.94 0.94 0.68 17 0.0087 4 0.83 0.68
##
## lower alpha upper 95% confidence boundaries
## 0.93 0.94 0.96
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## FFMQ2.3 0.93 0.94 0.93 0.67 15 0.0102 0.0039 0.67
## FFMQ2.10 0.94 0.94 0.93 0.68 15 0.0099 0.0026 0.69
## FFMQ2.14 0.94 0.94 0.94 0.70 16 0.0094 0.0021 0.69
## FFMQ2.17 0.94 0.94 0.93 0.68 15 0.0098 0.0031 0.69
## FFMQ2.25 0.93 0.93 0.93 0.66 14 0.0108 0.0019 0.66
## FFMQ2.30 0.94 0.94 0.94 0.68 15 0.0101 0.0038 0.68
## FFMQ2.35 0.93 0.94 0.93 0.67 14 0.0103 0.0039 0.66
## FFMQ2.39 0.94 0.94 0.94 0.68 15 0.0098 0.0030 0.68
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## FFMQ2.3 94 0.86 0.86 0.84 0.81 3.9 1.00
## FFMQ2.10 94 0.84 0.84 0.82 0.79 4.0 0.98
## FFMQ2.14 94 0.78 0.79 0.75 0.72 4.2 0.87
## FFMQ2.17 94 0.84 0.84 0.81 0.78 3.7 1.05
## FFMQ2.25 94 0.90 0.90 0.90 0.87 3.9 0.97
## FFMQ2.30 94 0.85 0.85 0.83 0.80 4.1 0.97
## FFMQ2.35 94 0.86 0.87 0.85 0.82 4.1 0.93
## FFMQ2.39 94 0.84 0.83 0.81 0.78 4.0 1.07
##
## Non missing response frequency for each item
## 1 2 3 3.94736842105263 4 5 miss
## FFMQ2.3 0.01 0.06 0.31 0.00 0.28 0.34 0
## FFMQ2.10 0.01 0.05 0.26 0.00 0.30 0.38 0
## FFMQ2.14 0.00 0.03 0.19 0.00 0.31 0.47 0
## FFMQ2.17 0.01 0.14 0.29 0.00 0.30 0.27 0
## FFMQ2.25 0.01 0.05 0.28 0.00 0.31 0.35 0
## FFMQ2.30 0.01 0.04 0.21 0.00 0.27 0.47 0
## FFMQ2.35 0.00 0.06 0.20 0.01 0.32 0.40 0
## FFMQ2.39 0.01 0.09 0.24 0.00 0.21 0.45 0psych::alpha(m2[2:40][non_react])##
## Reliability analysis
## Call: psych::alpha(x = m2[2:40][non_react])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.85 0.85 0.85 0.44 5.6 0.024 3.3 0.55 0.4
##
## lower alpha upper 95% confidence boundaries
## 0.8 0.85 0.89
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## FFMQ2.4 0.84 0.84 0.84 0.47 5.2 0.025 0.0124 0.48
## FFMQ2.9 0.82 0.82 0.82 0.43 4.6 0.028 0.0146 0.38
## FFMQ2.19 0.81 0.81 0.81 0.42 4.3 0.030 0.0116 0.38
## FFMQ2.21 0.85 0.85 0.84 0.48 5.5 0.025 0.0119 0.50
## FFMQ2.24 0.83 0.83 0.82 0.45 4.9 0.027 0.0102 0.41
## FFMQ2.29 0.82 0.82 0.81 0.43 4.6 0.029 0.0116 0.40
## FFMQ2.33 0.81 0.81 0.80 0.42 4.4 0.030 0.0092 0.38
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## FFMQ2.4 94 0.64 0.65 0.56 0.51 3.1 0.73
## FFMQ2.9 94 0.75 0.76 0.70 0.65 3.3 0.75
## FFMQ2.19 94 0.80 0.79 0.76 0.70 3.4 0.77
## FFMQ2.21 94 0.61 0.62 0.51 0.47 3.5 0.71
## FFMQ2.24 94 0.71 0.70 0.65 0.58 3.4 0.83
## FFMQ2.29 94 0.75 0.75 0.71 0.65 3.2 0.73
## FFMQ2.33 94 0.79 0.79 0.77 0.70 3.2 0.78
##
## Non missing response frequency for each item
## 1 2 3 4 4.05263157894737 4.73684210526316 5 miss
## FFMQ2.4 0.03 0.10 0.62 0.22 0.01 0.00 0.02 0
## FFMQ2.9 0.01 0.10 0.50 0.35 0.00 0.00 0.04 0
## FFMQ2.19 0.01 0.09 0.46 0.39 0.00 0.00 0.05 0
## FFMQ2.21 0.01 0.02 0.50 0.39 0.00 0.01 0.06 0
## FFMQ2.24 0.03 0.06 0.48 0.36 0.00 0.00 0.06 0
## FFMQ2.29 0.01 0.12 0.56 0.28 0.00 0.00 0.03 0
## FFMQ2.33 0.02 0.13 0.52 0.30 0.00 0.00 0.03 0psych::alpha(m2[2:40][observe])##
## Reliability analysis
## Call: psych::alpha(x = m2[2:40][observe])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.89 0.9 0.9 0.52 8.7 0.016 3.8 0.73 0.51
##
## lower alpha upper 95% confidence boundaries
## 0.86 0.89 0.93
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## FFMQ2.1 0.89 0.89 0.90 0.53 8.0 0.018 0.016 0.51
## FFMQ2.6 0.88 0.88 0.89 0.51 7.3 0.020 0.016 0.51
## FFMQ2.11 0.89 0.90 0.90 0.55 8.6 0.017 0.014 0.57
## FFMQ2.15 0.87 0.87 0.87 0.49 6.7 0.021 0.012 0.50
## FFMQ2.20 0.87 0.87 0.88 0.50 6.9 0.020 0.011 0.50
## FFMQ2.26 0.88 0.88 0.89 0.52 7.7 0.018 0.014 0.51
## FFMQ2.31 0.88 0.88 0.89 0.52 7.5 0.019 0.014 0.51
## FFMQ2.36 0.89 0.89 0.89 0.54 8.1 0.018 0.016 0.52
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## FFMQ2.1 94 0.71 0.71 0.65 0.61 3.6 0.96
## FFMQ2.6 94 0.80 0.80 0.76 0.72 3.5 1.08
## FFMQ2.11 94 0.65 0.65 0.58 0.54 3.6 0.93
## FFMQ2.15 94 0.87 0.87 0.87 0.82 4.1 0.86
## FFMQ2.20 94 0.84 0.84 0.84 0.79 3.9 0.91
## FFMQ2.26 94 0.75 0.75 0.71 0.66 4.1 1.02
## FFMQ2.31 94 0.77 0.76 0.73 0.68 4.2 0.99
## FFMQ2.36 94 0.70 0.71 0.65 0.60 3.7 0.91
##
## Non missing response frequency for each item
## 1 2 3 3.21052631578947 3.40540540540541 4 5 miss
## FFMQ2.1 0.02 0.11 0.33 0.00 0.00 0.38 0.16 0
## FFMQ2.6 0.04 0.14 0.26 0.00 0.01 0.36 0.19 0
## FFMQ2.11 0.02 0.11 0.29 0.00 0.00 0.45 0.14 0
## FFMQ2.15 0.01 0.03 0.16 0.00 0.00 0.45 0.35 0
## FFMQ2.20 0.01 0.05 0.22 0.00 0.00 0.41 0.30 0
## FFMQ2.26 0.01 0.09 0.16 0.00 0.00 0.31 0.44 0
## FFMQ2.31 0.03 0.02 0.14 0.01 0.00 0.31 0.49 0
## FFMQ2.36 0.02 0.04 0.36 0.00 0.00 0.38 0.19 0Model 1: v ~ valence * session
d1m is the MAP est df created above - rename d1me
d1me <- d1m–Valence: Baseline effects–
Depression
Put MAP ests, DASS scores, and endorsement differences at baseline into the same dataframe
# Match the IDs in the qairre and map ests datasets
dd1_red <- dd1 %>%
filter(ID %in% intersect(d1me$ID, dd1$ID)) %>% arrange(ID)
d1_red <- d1me %>% filter(ID %in% intersect(d1me$ID, dd1_red$ID)) %>% arrange(ID)# Combine them
if (all(d1_red$ID == dd1_red$ID)) fdf_depr <- data.frame(d1_red, dd1_red$dass_sum_score)Recalc baseline positive and negative endorsement diffs
nends_b <- bx_df %>% filter(session=="Pre", valence=="negative") %>%
group_by(subj_idx) %>% summarize(m=mean(response), mrt=mean(rt))
pends_b <- bx_df %>% filter(session=="Pre", valence=="positive") %>% group_by(subj_idx) %>%
summarize(m=mean(response), mrt=mean(rt))
if (all(pends_b$subj_idx==nends_b$subj_idx)) {
pends_b$pn_diff <- pends_b$m-nends_b$m
pends_b$pn_rt_diff <- pends_b$mrt-nends_b$mrt
}# Match and combine the endorsements dataset with the depression one
ends_sub <- pends_b %>% filter(subj_idx %in% fdf_depr$ID)
# Add the difference in positive-negative endorsements to the depression dataset
if (all(ends_sub$subj_idx==fdf_depr$ID)) fdf_depr <- data.frame(fdf_depr, "pos_neg_ends_diff"=ends_sub$pn_diff)rt <- bx_df %>%
group_by(subj_idx) %>% summarize(m=mean(rt), mrt=mean(rt))
rt_val <- bx_df %>%
group_by(subj_idx, valence) %>% summarize(mrt=mean(rt))## `summarise()` has grouped output by 'subj_idx'. You can override using the
## `.groups` argument.rt_val_short <- data.frame(
"ID"=unique(rt_val$subj_idx),
"overall_rt_diff"=rt_val[rt_val$valence == "positive", "mrt"]-rt_val[rt_val$valence == "negative", "mrt"])
nends_p <- bx_df %>% filter(session=="Post", valence=="negative") %>%
group_by(subj_idx) %>% summarize(m=mean(response), mrt=mean(rt))
pends_p <- bx_df %>% filter(session=="Post", valence=="positive") %>% group_by(subj_idx) %>% summarize(m=mean(response), mrt=mean(rt))
if (all(pends_p$subj_idx==nends_p$subj_idx)) {
pends_p$pn_diff <- pends_p$m-nends_p$m
pends_p$pn_rt_diff <- pends_p$mrt-nends_p$mrt
}Depression relationship
cor.test(fdf_depr$dd1_red.dass_sum_score, fdf_depr$v_val_ctr_)##
## Pearson's product-moment correlation
##
## data: fdf_depr$dd1_red.dass_sum_score and fdf_depr$v_val_ctr_
## t = -4.0199, df = 93, p-value = 0.0001182
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5441034 -0.1986105
## sample estimates:
## cor
## -0.3847524cor.test(fdf_depr$dd1_red.dass_sum_score, fdf_depr$pos_neg_ends_diff)##
## Pearson's product-moment correlation
##
## data: fdf_depr$dd1_red.dass_sum_score and fdf_depr$pos_neg_ends_diff
## t = -4.4162, df = 93, p-value = 2.709e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5700670 -0.2344955
## sample estimates:
## cor
## -0.4163609Figure 4B, right
dvv <- ggplot(fdf_depr, aes(v_val_ctr_, dd1_red.dass_sum_score)) +
geom_point(size=8, alpha=1, pch=21, fill="gray57") +
#ggtitle(paste0( "r=-.30 p < .001 ")) +
geom_smooth(method="lm", color="black") +
ga + ap +
geom_smooth(method="lm", color="black") +
ga + theme(plot.title = element_text(margin=margin(b=-3),
size=30, hjust=.15, vjust=.6)) + theme(axis.title = element_text(size=35), axis.text=element_text(size=35)) +
ylab("baseline depression \n (via DASS)") + xlab("positive > negative \n drift rate estimates") +
theme(plot.title = element_text(margin=margin(b=5),
size=27, hjust=.15, vjust=.6)) +
geom_text(x=2.1, y=24.5, size=15, label="*****") + ylim(0, 37)
dvv## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
#ggsave("../paper/figs/depr_v-val.png", dvv, width=9, height=7)Mindfulness
# Match the IDs in the qairre and map ests datasets
m1_r <- m1 %>%
filter(ID %in% intersect(d1me$ID, m1$ID)) %>% arrange(ID)
d1_r <- d1me %>% filter(ID %in% intersect(d1me$ID, m1_r$ID)) %>% arrange(ID)
# Combine them
if (all(d1_r$ID == m1_r$ID)) fdf <- data.frame(d1_r,
m1_r$aware_sum, m1_r$nonjudge_sum, m1_r$nonreact_sum, m1_r$observe_sum)Add the endorsement difference variable
ends_sub_m <- pends_b %>% filter(subj_idx %in% m1_r$ID)
if (all(m1_r$ID == ends_sub_m$subj_idx)) fdf <- data.frame(fdf, "pos_neg_diff_m"=ends_sub_m$pn_diff)summary(
md_ddm_1 <- lm(v_val_ctr_ ~
scale(m1_r.aware_sum)*scale(m1_r.nonjudge_sum) +
scale(m1_r.nonreact_sum) ,
data=fdf)
)##
## Call:
## lm(formula = v_val_ctr_ ~ scale(m1_r.aware_sum) * scale(m1_r.nonjudge_sum) +
## scale(m1_r.nonreact_sum), data = fdf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70168 -0.26863 -0.02867 0.23972 0.91191
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.322348 0.039716 33.295
## scale(m1_r.aware_sum) 0.104288 0.041960 2.485
## scale(m1_r.nonjudge_sum) 0.050452 0.045052 1.120
## scale(m1_r.nonreact_sum) 0.019338 0.044843 0.431
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum) 0.004105 0.036826 0.111
## Pr(>|t|)
## (Intercept) <2e-16 ***
## scale(m1_r.aware_sum) 0.0148 *
## scale(m1_r.nonjudge_sum) 0.2658
## scale(m1_r.nonreact_sum) 0.6673
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum) 0.9115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.359 on 90 degrees of freedom
## Multiple R-squared: 0.1432, Adjusted R-squared: 0.1052
## F-statistic: 3.761 on 4 and 90 DF, p-value: 0.007092car::vif(md_ddm_1)## scale(m1_r.aware_sum)
## 1.284368
## scale(m1_r.nonjudge_sum)
## 1.480637
## scale(m1_r.nonreact_sum)
## 1.466887
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum)
## 1.164014Compare to raw endorsement difference as outcome var
summary(
md_raw_1 <- lm(pos_neg_diff_m ~
scale(m1_r.aware_sum)*scale(m1_r.nonjudge_sum) +
scale(m1_r.nonreact_sum) ,
data=fdf)
)##
## Call:
## lm(formula = pos_neg_diff_m ~ scale(m1_r.aware_sum) * scale(m1_r.nonjudge_sum) +
## scale(m1_r.nonreact_sum), data = fdf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.54962 -0.08583 -0.00330 0.14273 0.32519
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.714502 0.019190 37.233
## scale(m1_r.aware_sum) 0.063934 0.020275 3.153
## scale(m1_r.nonjudge_sum) 0.023046 0.021769 1.059
## scale(m1_r.nonreact_sum) 0.001903 0.021667 0.088
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum) -0.008999 0.017794 -0.506
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## scale(m1_r.aware_sum) 0.00219 **
## scale(m1_r.nonjudge_sum) 0.29259
## scale(m1_r.nonreact_sum) 0.93020
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum) 0.61429
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1734 on 90 degrees of freedom
## Multiple R-squared: 0.1684, Adjusted R-squared: 0.1314
## F-statistic: 4.556 on 4 and 90 DF, p-value: 0.002135car::vif(md_raw_1)## scale(m1_r.aware_sum)
## 1.284368
## scale(m1_r.nonjudge_sum)
## 1.480637
## scale(m1_r.nonreact_sum)
## 1.466887
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum)
## 1.164014out <- summary(md_ddm_1)
reg_df <- data.table(c("aware", "non-judge", "non-react", "aware*non-judge"),
out$coefficients[2:5, 1],
out$coefficients[2:5, 2]) %>% setNames(c("coefficient", "estimate", "se"))
reg_df$coefficient <- factor(reg_df$coefficient, levels=c("aware", "non-judge", "non-react", "aware*non-judge"))mr <- ggplot(reg_df, aes(x=coefficient, y=estimate)) +
geom_errorbar(aes(ymin=estimate-se, ymax=estimate+se), size=2, width=0) +
geom_point(fill="gray57", color="black", alpha=.8, size=8, pch=21) +
geom_hline(yintercept=c(0), linetype="dotted") +
theme(axis.text.x = element_text(angle=25, hjust=1)) +
ga + ap + ylab("coefficients (±1 SEM)") + xlab("") +
ggtitle("Regression coefficients \n of pos > neg drift rate ~ FFMQ mindfulness facets") +
theme(plot.title = element_text(margin=margin(b=5),
size=27, hjust=.15, vjust=.6)) +
geom_text(x=1, y=.17, size=15, label="*") + ylim(-.05, .18) +
theme(axis.text = element_text(size=30),
axis.title = element_text(size=30))
mr
#ggsave("../paper/figs/mindful_regr.png", mr, width=13, height=6)Adding observe
summary(
md_wobs_1 <- lm(v_val_ctr_ ~
scale(m1_r.aware_sum)*scale(m1_r.nonjudge_sum) +
scale(m1_r.nonreact_sum) + scale(m1_r.observe_sum),
data=fdf)
)##
## Call:
## lm(formula = v_val_ctr_ ~ scale(m1_r.aware_sum) * scale(m1_r.nonjudge_sum) +
## scale(m1_r.nonreact_sum) + scale(m1_r.observe_sum), data = fdf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70167 -0.26870 -0.02481 0.24000 0.90945
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.322071 0.040042 33.017
## scale(m1_r.aware_sum) 0.104159 0.042215 2.467
## scale(m1_r.nonjudge_sum) 0.051018 0.045690 1.117
## scale(m1_r.nonreact_sum) 0.018006 0.047208 0.381
## scale(m1_r.observe_sum) 0.003742 0.039247 0.095
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum) 0.004791 0.037723 0.127
## Pr(>|t|)
## (Intercept) <2e-16 ***
## scale(m1_r.aware_sum) 0.0155 *
## scale(m1_r.nonjudge_sum) 0.2672
## scale(m1_r.nonreact_sum) 0.7038
## scale(m1_r.observe_sum) 0.9243
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum) 0.8992
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.361 on 89 degrees of freedom
## Multiple R-squared: 0.1433, Adjusted R-squared: 0.09519
## F-statistic: 2.978 on 5 and 89 DF, p-value: 0.01564car::vif(md_wobs_1)## scale(m1_r.aware_sum)
## 1.285691
## scale(m1_r.nonjudge_sum)
## 1.506095
## scale(m1_r.nonreact_sum)
## 1.607810
## scale(m1_r.observe_sum)
## 1.111253
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum)
## 1.207987summary(
md_wobs_2 <- lm(v_val_ctr_ ~
scale(m1_r.aware_sum)*scale(m1_r.nonjudge_sum) +
scale(m1_r.nonreact_sum) + scale(m1_r.observe_sum)*scale(m1_r.nonjudge_sum),
data=fdf)
)##
## Call:
## lm(formula = v_val_ctr_ ~ scale(m1_r.aware_sum) * scale(m1_r.nonjudge_sum) +
## scale(m1_r.nonreact_sum) + scale(m1_r.observe_sum) * scale(m1_r.nonjudge_sum),
## data = fdf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.69522 -0.27415 -0.03104 0.23719 0.89794
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.319376 0.040230 32.796
## scale(m1_r.aware_sum) 0.108939 0.042654 2.554
## scale(m1_r.nonjudge_sum) 0.055135 0.046018 1.198
## scale(m1_r.nonreact_sum) 0.006982 0.049032 0.142
## scale(m1_r.observe_sum) -0.004312 0.040437 -0.107
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum) 0.008756 0.038070 0.230
## scale(m1_r.nonjudge_sum):scale(m1_r.observe_sum) 0.030466 0.035884 0.849
## Pr(>|t|)
## (Intercept) <2e-16 ***
## scale(m1_r.aware_sum) 0.0124 *
## scale(m1_r.nonjudge_sum) 0.2341
## scale(m1_r.nonreact_sum) 0.8871
## scale(m1_r.observe_sum) 0.9153
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum) 0.8186
## scale(m1_r.nonjudge_sum):scale(m1_r.observe_sum) 0.3982
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3615 on 88 degrees of freedom
## Multiple R-squared: 0.1503, Adjusted R-squared: 0.09234
## F-statistic: 2.594 on 6 and 88 DF, p-value: 0.02317car::vif(md_wobs_2)## scale(m1_r.aware_sum)
## 1.308487
## scale(m1_r.nonjudge_sum)
## 1.523001
## scale(m1_r.nonreact_sum)
## 1.729043
## scale(m1_r.observe_sum)
## 1.175969
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum)
## 1.226443
## scale(m1_r.nonjudge_sum):scale(m1_r.observe_sum)
## 1.202048Robustness check: Confirm effects are same with casewise deletion
# Match the IDs in the qairre and map ests datasets
m1_r_cd <- m1_casewise_delete %>%
filter(ID %in% intersect(d1me$ID, m1_casewise_delete$ID)) %>% arrange(ID)
d1_r_cd <- d1me %>% filter(ID %in% intersect(d1me$ID, m1_casewise_delete$ID)) %>% arrange(ID)
# Combine them
if (all(d1_r_cd$ID == m1_r_cd$ID)) fdf_cd <- data.frame(d1_r_cd,
m1_r_cd$aware_sum,
m1_r_cd$nonreact_sum,
m1_r_cd$nonjudge_sum)summary(
md_d_1_cd <- lm(v_val_ctr_ ~
scale(m1_r_cd.aware_sum)*scale(m1_r_cd.nonjudge_sum) +
scale(m1_r_cd.nonreact_sum),
data=fdf_cd)
)##
## Call:
## lm(formula = v_val_ctr_ ~ scale(m1_r_cd.aware_sum) * scale(m1_r_cd.nonjudge_sum) +
## scale(m1_r_cd.nonreact_sum), data = fdf_cd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.74755 -0.24873 -0.02561 0.22460 0.87001
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 1.365236 0.040250
## scale(m1_r_cd.aware_sum) 0.101092 0.042752
## scale(m1_r_cd.nonjudge_sum) 0.045756 0.045607
## scale(m1_r_cd.nonreact_sum) 0.037034 0.046448
## scale(m1_r_cd.aware_sum):scale(m1_r_cd.nonjudge_sum) -0.006047 0.037442
## t value Pr(>|t|)
## (Intercept) 33.919 <2e-16 ***
## scale(m1_r_cd.aware_sum) 2.365 0.0204 *
## scale(m1_r_cd.nonjudge_sum) 1.003 0.3187
## scale(m1_r_cd.nonreact_sum) 0.797 0.4276
## scale(m1_r_cd.aware_sum):scale(m1_r_cd.nonjudge_sum) -0.162 0.8721
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3487 on 81 degrees of freedom
## Multiple R-squared: 0.1557, Adjusted R-squared: 0.114
## F-statistic: 3.735 on 4 and 81 DF, p-value: 0.007702car::vif(md_d_1_cd)## scale(m1_r_cd.aware_sum)
## 1.277384
## scale(m1_r_cd.nonjudge_sum)
## 1.453664
## scale(m1_r_cd.nonreact_sum)
## 1.507804
## scale(m1_r_cd.aware_sum):scale(m1_r_cd.nonjudge_sum)
## 1.170969fd_r <- fdf %>% filter(ID %in% dd1_casewise_delete$ID) #%>% select(v_val_ctr_)
d1_r <- dd1_casewise_delete %>% filter(ID %in% fd_r$ID)
if (all(d1_r$ID==fd_r$ID)) cor.test(fd_r$v_val_ctr_, d1_r$dass_sum_score)##
## Pearson's product-moment correlation
##
## data: fd_r$v_val_ctr_ and d1_r$dass_sum_score
## t = -3.387, df = 90, p-value = 0.00105
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5061892 -0.1411328
## sample estimates:
## cor
## -0.336231Mindfulness and depression in same
if (all(fdf$ID == fdf_depr$ID)) fdf <- data.frame(fdf, fdf_depr$dd1_red.dass_sum_score)Depression relationship
# Sanity check
#cor.test(fdf$fdf_depr.dd1_red.dass_sum_score, fdf$v_val_ctr_)
cor.test(fdf$fdf_depr.dd1_red.dass_sum_score, fdf$m1_r.aware_sum)##
## Pearson's product-moment correlation
##
## data: fdf$fdf_depr.dd1_red.dass_sum_score and fdf$m1_r.aware_sum
## t = -4.4046, df = 93, p-value = 2.832e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5693236 -0.2334552
## sample estimates:
## cor
## -0.4154507Check other mindfulness facet - depression relationships
- Comparably strong relationship with non-judgment but weaker with non-react
cor.test(fdf$fdf_depr.dd1_red.dass_sum_score, fdf$m1_r.nonjudge_sum)##
## Pearson's product-moment correlation
##
## data: fdf$fdf_depr.dd1_red.dass_sum_score and fdf$m1_r.nonjudge_sum
## t = -4.0497, df = 93, p-value = 0.0001061
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5461076 -0.2013479
## sample estimates:
## cor
## -0.3871791cor.test(fdf$fdf_depr.dd1_red.dass_sum_score, fdf$m1_r.nonreact_sum)##
## Pearson's product-moment correlation
##
## data: fdf$fdf_depr.dd1_red.dass_sum_score and fdf$m1_r.nonreact_sum
## t = -1.4622, df = 93, p-value = 0.1471
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3411428 0.0532451
## sample estimates:
## cor
## -0.1499066summary(
all_m <- lm(v_val_ctr_ ~
scale(m1_r.aware_sum)*scale(m1_r.nonjudge_sum) +
scale(m1_r.nonreact_sum) + scale(fdf_depr.dd1_red.dass_sum_score),
data=fdf)
)##
## Call:
## lm(formula = v_val_ctr_ ~ scale(m1_r.aware_sum) * scale(m1_r.nonjudge_sum) +
## scale(m1_r.nonreact_sum) + scale(fdf_depr.dd1_red.dass_sum_score),
## data = fdf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7122 -0.2505 -0.0240 0.2466 0.8261
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.31119 0.03862 33.955
## scale(m1_r.aware_sum) 0.06121 0.04358 1.405
## scale(m1_r.nonjudge_sum) 0.02399 0.04464 0.537
## scale(m1_r.nonreact_sum) 0.02105 0.04335 0.485
## scale(fdf_depr.dd1_red.dass_sum_score) -0.11530 0.04266 -2.703
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum) 0.03175 0.03704 0.857
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## scale(m1_r.aware_sum) 0.16362
## scale(m1_r.nonjudge_sum) 0.59239
## scale(m1_r.nonreact_sum) 0.62855
## scale(fdf_depr.dd1_red.dass_sum_score) 0.00823 **
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum) 0.39365
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.347 on 89 degrees of freedom
## Multiple R-squared: 0.2082, Adjusted R-squared: 0.1637
## F-statistic: 4.681 on 5 and 89 DF, p-value: 0.0007688car::vif(all_m)## scale(m1_r.aware_sum)
## 1.482614
## scale(m1_r.nonjudge_sum)
## 1.555482
## scale(m1_r.nonreact_sum)
## 1.467199
## scale(fdf_depr.dd1_red.dass_sum_score)
## 1.420407
## scale(m1_r.aware_sum):scale(m1_r.nonjudge_sum)
## 1.260104Univariate
cor.test(fdf$v_val_ctr_, fdf$m1_r.aware_sum)##
## Pearson's product-moment correlation
##
## data: fdf$v_val_ctr_ and fdf$m1_r.aware_sum
## t = 3.5931, df = 93, p-value = 0.0005249
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1587662 0.5144785
## sample estimates:
## cor
## 0.3491374Change over time regressions
Add dosage for time
# big_df <- haven::read_sav("./../../data/raw_files/K23_Hitchcock_medpracticedata.sav")
# big_df <- big_df %>% filter(ID %in% fdf$ID) %>% arrange(ID)
# fdf_md <- fdf %>% filter(ID %in% big_df$ID) %>% arrange(ID)
# if (all(fdf_md$ID==big_df$ID)) cor.test(fdf_md$v_val_ctr.sess_ctr_, big_df$TOTALminpwerwk.8wks)Recreate post - pre scores
# Post scores
nends_p <- bx_df %>% filter(session=="Post", valence=="negative") %>%
group_by(subj_idx) %>% summarize(m=mean(response), mrt=mean(rt))
pends_p <- bx_df %>% filter(session=="Post", valence=="positive") %>% group_by(subj_idx) %>% summarize(m=mean(response), mrt=mean(rt))
if (all(pends_p$subj_idx==nends_p$subj_idx)) {
pends_p$pn_diff <- pends_p$m-nends_p$m
}Create a difference variable of (post pos - neg) - (pre pos - neg)
if (all(pends_p$subj_idx==pends_b$subj_idx)) pends_p$diff_in_diff_pn_ends <- pends_p$pn_diff - pends_b$pn_diffWe are just looking at FFMQ and depression at pre- and post-intervention and want to relate these to the SRET valence*time regressor. Had originally tried mixed-effects models to try to get time REs, but had convergence issues which make sense given low data per pt, so opting for simpler approach of change scores
m1_shared <- m1 %>% filter(ID %in% m2$ID)
m2_shared <- m2 %>% filter(ID %in% m1$ID)
if (all(m1_shared$ID==m2_shared$ID)) {
mindfulness_changes <- data.frame(m1_shared$ID,
m2_shared$aware_sum-m1_shared$aware_sum,
m2_shared$nonjudge_sum-m1_shared$nonjudge_sum,
m2_shared$nonreact_sum-m1_shared$nonreact_sum,
m2_shared$observe_sum-m1_shared$observe_sum) %>%
setNames(c("ID", "aware_change", "nj_change", "nr_change", "obs_change"))
}
hist(mindfulness_changes$nr_change)
hist(mindfulness_changes$nj_change)
hist(mindfulness_changes$aware_change)
hist(mindfulness_changes$obs_change)
#hist(fdf$FFMQ_Aware_SUM_Q1)
#hist(fdf$FFMQ_Nonjudge_SUM_Q1)
a <- ggplot(fdf, aes(x=m1_r.aware_sum)) +
geom_histogram(bins=25, fill="gray38", color="black") +
xlim(5, 45) + ga + theme(axis.text=element_blank(), axis.ticks = element_blank()) +
ylab("") + xlab("") + ggtitle("Distributions of Mindfulness Facets at Baseline") +
theme(plot.title = element_text(margin=margin(b=-5),
size=35, hjust=.15, vjust=.6)) +
annotate("text", x = 40, y = 10, label = "Aware", size=10)
# annotate("text", x = 2.45, y = .85, label = TeX("yellow: $\\alpha$ < median"), size=13)
b <- ggplot(fdf, aes(x=m1_r.nonjudge_sum)) +
geom_histogram(bins=25, fill="gray38", color="black") +
xlim(5, 45) + ga + ap + ylab("") + xlab("") +
theme(axis.text.y=element_blank(), axis.ticks.y = element_blank()) +
annotate("text", x = 35, y = 10, label = "Non-judgmental", size=10)
comb_m <- a/b
comb_m## Warning: Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
#ggsave("../paper/figs/mind_dist.png", comb_m, width=18, height=8)# mindfulness_changes$aware_change
# mindfulness_changes$aware_change
# range(mindfulness_changes$aware_change)
# range(mindfulness_changes$nj_change)
c <- ggplot(mindfulness_changes, aes(x=aware_change)) +
geom_histogram(bins=25, fill="gray75", color="black") +
ga + theme(axis.text=element_blank(), axis.ticks = element_blank()) +
ylab("") + xlab("") + ggtitle("Distributions of Change (Post - Pre) in Mindfulness Facets") +
theme(plot.title = element_text(margin=margin(b=-5),
size=35, hjust=.15, vjust=.6)) +
xlim(-14, 27) +
annotate("text", x = 17, y = 11, label = "Aware", size=10)
# annotate("text", x = 2.45, y = .85, label = TeX("yellow: $\\alpha$ < median"), size=13)
d <- ggplot(mindfulness_changes, aes(x=nj_change)) +
geom_histogram(bins=25, fill="gray75", color="black") +
ga + ap + ylab("") + xlab("") +
theme(axis.text.y=element_blank(), axis.ticks.y = element_blank()) +
xlim(-14, 27) +
annotate("text", x = 17, y = 15, label = "Non-judgmental", size=10)
comb_ch <- c/d
comb_ch## Warning: Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
#ggsave("../paper/figs/mind_ch.png", comb_ch, width=18, height=8)d1me_shared <- d1me %>% filter(ID %in% mindfulness_changes$ID)Examining how the machine learning procedure by Beevers et al. 19 JAP influences sparsity. Model draws on code included in the paper
var_names <- c(
# 3, 5, 10, 13
"no-positive-feelings", "couldnt-get-going", "easily-upset", "sad-and-depressed",
# 16, 17, 21
"lost-interest", "not-worth-much-as-person", "life-not-worthwhile",
# 24, 26, 31, 34, 37,
"no-enjoyment", "down-hearted-and-blue", "no-enthusiasm", "worthless", "no-hope",
#38, 42
"life-meaningless", "no-initiative",
"drift_regressor"
)Empirical
if (all(d1_red$ID == dd1_red$ID)) {
t1_depr_combined <- data.frame(dd1_red, "v_val_ctr"=d1_red$v_val_ctr_)
}
#cor.test(t1_depr_combined$v_val_ctr, t1_depr_combined$dass_sum_score)
cmat <- cor(t1_depr_combined %>% select(-c(ID, dass_sum_score)))
#cmatjust_rel_vars_emp <- cbind(t1_depr_combined[grep("DASS1", names(t1_depr_combined))], t1_depr_combined$v_val_ctr) %>%
setNames(var_names)
rf_out_emp <- beset_rf(drift_regressor ~ ., data=just_rel_vars_emp)
importance(rf_out_emp)
To examine the influence of sample size
(commented because time consuming)
# sample_sizes <- c(100, 200, 500, 1e3, 5e3)
#
#
# for (ss in seq_along(sample_sizes)) {
# simulated_data <- rnorm_multi(n=sample_sizes[ss], vars=ncol(cmat), mu=0, sd=1, r=cmat)
# # Impose constraints on the depression items
# depr_items <- round(simulated_data[1:14, 1:14], 0)
# depr_items[depr_items < 0] <- 0
# depr_items[depr_items > 3] <- 3
#
# # Put them back in
# simulated_data[1:14, 1:14] <- depr_items
# simulated_data <- simulated_data %>% setNames(var_names)
#
# rf_out <- beset_rf(drift_regressor ~ ., data=simulated_data)
# print(importance(rf_out))
# }Mindfulness change
# Add drift rate change regressor
if (all(d1me_shared$ID==mindfulness_changes$ID)) {
mindfulness_changes$dr_change <- d1me_shared$v_val_ctr.sess_ctr_
}Add scaled variables
mindfulness_changes <-
data.frame(mindfulness_changes,
data.frame(
"aware_ch_sc"=scale(mindfulness_changes$aware_change),
"nj_ch_sc"=scale(mindfulness_changes$nj_change),
"nr_ch_sc"=scale(mindfulness_changes$nr_change))
)Add the difference in difference endorsements variable to the df with mindfulness changes
diff_ends_sub <- pends_p %>% filter(subj_idx %in% mindfulness_changes$ID)
if (all(diff_ends_sub$subj_idx==mindfulness_changes$ID)) mindfulness_changes<- data.frame(mindfulness_changes, "diff_in_diff_posneg"=diff_ends_sub$diff_in_diff_pn_ends)summary(
change_1m <-
lm(dr_change ~
aware_ch_sc*nj_ch_sc +
nr_ch_sc,
data=mindfulness_changes)
)##
## Call:
## lm(formula = dr_change ~ aware_ch_sc * nj_ch_sc + nr_ch_sc, data = mindfulness_changes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50841 -0.11252 -0.01161 0.11116 0.46922
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.128739 0.020388 6.315 1.07e-08 ***
## aware_ch_sc 0.060368 0.022153 2.725 0.00776 **
## nj_ch_sc 0.025647 0.023174 1.107 0.27143
## nr_ch_sc -0.030873 0.021720 -1.421 0.15872
## aware_ch_sc:nj_ch_sc -0.007739 0.016166 -0.479 0.63332
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1808 on 88 degrees of freedom
## Multiple R-squared: 0.133, Adjusted R-squared: 0.09361
## F-statistic: 3.375 on 4 and 88 DF, p-value: 0.01283car::vif(change_1m)## aware_ch_sc nj_ch_sc nr_ch_sc
## 1.381155 1.511342 1.327653
## aware_ch_sc:nj_ch_sc
## 1.034201summary(
change_raw_m <-
lm(diff_in_diff_posneg ~
aware_ch_sc*nj_ch_sc +
nr_ch_sc,
data=mindfulness_changes)
)##
## Call:
## lm(formula = diff_in_diff_posneg ~ aware_ch_sc * nj_ch_sc + nr_ch_sc,
## data = mindfulness_changes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53084 -0.08525 -0.00655 0.09687 0.47572
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.070107 0.017923 3.912 0.00018 ***
## aware_ch_sc 0.043281 0.019475 2.222 0.02882 *
## nj_ch_sc 0.008860 0.020372 0.435 0.66470
## nr_ch_sc -0.003591 0.019094 -0.188 0.85124
## aware_ch_sc:nj_ch_sc 0.007342 0.014212 0.517 0.60672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1589 on 88 degrees of freedom
## Multiple R-squared: 0.08538, Adjusted R-squared: 0.04381
## F-statistic: 2.054 on 4 and 88 DF, p-value: 0.09373car::vif(change_raw_m)## aware_ch_sc nj_ch_sc nr_ch_sc
## 1.381155 1.511342 1.327653
## aware_ch_sc:nj_ch_sc
## 1.034201out_2 <- summary(change_1m)
out_2##
## Call:
## lm(formula = dr_change ~ aware_ch_sc * nj_ch_sc + nr_ch_sc, data = mindfulness_changes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50841 -0.11252 -0.01161 0.11116 0.46922
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.128739 0.020388 6.315 1.07e-08 ***
## aware_ch_sc 0.060368 0.022153 2.725 0.00776 **
## nj_ch_sc 0.025647 0.023174 1.107 0.27143
## nr_ch_sc -0.030873 0.021720 -1.421 0.15872
## aware_ch_sc:nj_ch_sc -0.007739 0.016166 -0.479 0.63332
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1808 on 88 degrees of freedom
## Multiple R-squared: 0.133, Adjusted R-squared: 0.09361
## F-statistic: 3.375 on 4 and 88 DF, p-value: 0.01283reg_df_2 <- data.table(c("aware-change", "non-judge-change", "non-react-change", "aware*non-judge-change"),
out_2$coefficients[2:5, 1],
out_2$coefficients[2:5, 2]) %>% setNames(c("coefficient", "estimate", "se"))
reg_df_2$coefficient <- factor(reg_df_2$coefficient, levels=c("aware-change", "non-judge-change", "non-react-change", "aware*non-judge-change"))mr_2 <- ggplot(reg_df_2, aes(x=coefficient, y=estimate)) +
geom_errorbar(aes(ymin=estimate-se, ymax=estimate+se), color="gray57", size=2, width=0) +
geom_point(fill="black", color="black", alpha=.3, size=8, pch=21) +
geom_hline(yintercept=c(0), linetype="dotted") +
theme(axis.text.x = element_text(angle=25, hjust=1)) +
ga + ap + ylab("coefficients (±1 SEM)") + xlab("") +
ggtitle(paste0("Regression coefficients \n of pos > neg drift rate * time ~ FFMQ mindfulness facets")) +
theme(plot.title = element_text(margin=margin(b=5),
size=27, hjust=.15, vjust=.6)) +
geom_text(x=1, y=.092, size=15, label="**") + ylim(-.08, .10) +
theme(axis.text = element_text(size=30),
axis.title = element_text(size=30)) +
scale_x_discrete(labels=c("aware-change" = "aware", "non-judge-change" = "non-judge",
"non-react-change" = "non-react", "aware*non-judge-change" = "aware*non-judge"))
mr_2
#ggsave("../paper/figs/mindful_regr_change.png", mr_2, width=13, height=6)Robustness check - add treatment condition
assign_short <- extra_df_full %>% filter(ID %in% mindfulness_changes$ID)testit::assert(assign_short$ID == mindfulness_changes$ID)
mindfulness_changes$tx_cond <- assign_short$Treatmentmindfulness_changes$tx_cond_f <- factor(mindfulness_changes$tx_cond, levels=c("OM", "FA", "MBCT"))
contrasts(mindfulness_changes$tx_cond_f) <- contr.sum(3)
#mindfulness_changes$tx_cond_fsummary(
change_1m_tx <-
lm(dr_change ~ tx_cond_f/(
scale(aware_change)*scale(nj_change) +
scale(nr_change)),
data=mindfulness_changes)) ##
## Call:
## lm(formula = dr_change ~ tx_cond_f/(scale(aware_change) * scale(nj_change) +
## scale(nr_change)), data = mindfulness_changes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.44785 -0.08492 -0.00807 0.10600 0.42981
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.121422 0.019360 6.272
## tx_cond_f1 -0.044635 0.027901 -1.600
## tx_cond_f2 -0.001769 0.026902 -0.066
## tx_cond_fOM:scale(aware_change) 0.108505 0.040931 2.651
## tx_cond_fFA:scale(aware_change) -0.009815 0.036175 -0.271
## tx_cond_fMBCT:scale(aware_change) 0.134140 0.034059 3.938
## tx_cond_fOM:scale(nj_change) -0.002977 0.041566 -0.072
## tx_cond_fFA:scale(nj_change) 0.052491 0.035741 1.469
## tx_cond_fMBCT:scale(nj_change) 0.002369 0.043907 0.054
## tx_cond_fOM:scale(nr_change) -0.033392 0.043222 -0.773
## tx_cond_fFA:scale(nr_change) -0.005319 0.026502 -0.201
## tx_cond_fMBCT:scale(nr_change) -0.018680 0.054729 -0.341
## tx_cond_fOM:scale(aware_change):scale(nj_change) -0.055346 0.024277 -2.280
## tx_cond_fFA:scale(aware_change):scale(nj_change) 0.028655 0.033664 0.851
## tx_cond_fMBCT:scale(aware_change):scale(nj_change) 0.045439 0.029203 1.556
## Pr(>|t|)
## (Intercept) 1.84e-08 ***
## tx_cond_f1 0.113688
## tx_cond_f2 0.947753
## tx_cond_fOM:scale(aware_change) 0.009718 **
## tx_cond_fFA:scale(aware_change) 0.786868
## tx_cond_fMBCT:scale(aware_change) 0.000177 ***
## tx_cond_fOM:scale(nj_change) 0.943078
## tx_cond_fFA:scale(nj_change) 0.145954
## tx_cond_fMBCT:scale(nj_change) 0.957113
## tx_cond_fOM:scale(nr_change) 0.442115
## tx_cond_fFA:scale(nr_change) 0.841443
## tx_cond_fMBCT:scale(nr_change) 0.733780
## tx_cond_fOM:scale(aware_change):scale(nj_change) 0.025354 *
## tx_cond_fFA:scale(aware_change):scale(nj_change) 0.397266
## tx_cond_fMBCT:scale(aware_change):scale(nj_change) 0.123763
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1664 on 78 degrees of freedom
## Multiple R-squared: 0.3493, Adjusted R-squared: 0.2326
## F-statistic: 2.991 on 14 and 78 DF, p-value: 0.001055Robustness check: Add observe change
summary(
change_0m <-
lm(dr_change ~
scale(aware_change)*scale(nj_change) +
scale(nr_change) +
scale(obs_change),
data=mindfulness_changes)
)##
## Call:
## lm(formula = dr_change ~ scale(aware_change) * scale(nj_change) +
## scale(nr_change) + scale(obs_change), data = mindfulness_changes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.51231 -0.11970 -0.00832 0.11091 0.47367
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.128903 0.020503 6.287 1.25e-08 ***
## scale(aware_change) 0.059524 0.022465 2.650 0.00957 **
## scale(nj_change) 0.026254 0.023393 1.122 0.26482
## scale(nr_change) -0.034735 0.025688 -1.352 0.17981
## scale(obs_change) 0.006682 0.023414 0.285 0.77603
## scale(aware_change):scale(nj_change) -0.008070 0.016293 -0.495 0.62162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1818 on 87 degrees of freedom
## Multiple R-squared: 0.1338, Adjusted R-squared: 0.08405
## F-statistic: 2.688 on 5 and 87 DF, p-value: 0.02624car::vif(change_0m)## scale(aware_change) scale(nj_change)
## 1.405536 1.523951
## scale(nr_change) scale(obs_change)
## 1.837643 1.526737
## scale(aware_change):scale(nj_change)
## 1.039469summary(
change_1m <-
lm(dr_change ~
scale(aware_change)*scale(nj_change) +
scale(nr_change) +
scale(obs_change)*scale(nj_change),
data=mindfulness_changes)
)##
## Call:
## lm(formula = dr_change ~ scale(aware_change) * scale(nj_change) +
## scale(nr_change) + scale(obs_change) * scale(nj_change),
## data = mindfulness_changes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.51286 -0.11085 -0.00314 0.11083 0.47533
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.130476 0.020740 6.291 1.27e-08 ***
## scale(aware_change) 0.059942 0.022558 2.657 0.00939 **
## scale(nj_change) 0.027343 0.023546 1.161 0.24875
## scale(nr_change) -0.039025 0.026730 -1.460 0.14795
## scale(obs_change) 0.007449 0.023533 0.317 0.75237
## scale(aware_change):scale(nj_change) -0.006046 0.016688 -0.362 0.71801
## scale(nj_change):scale(obs_change) -0.011261 0.018533 -0.608 0.54505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1824 on 86 degrees of freedom
## Multiple R-squared: 0.1375, Adjusted R-squared: 0.07736
## F-statistic: 2.286 on 6 and 86 DF, p-value: 0.0428car::vif(change_1m)## scale(aware_change) scale(nj_change)
## 1.406847 1.532839
## scale(nr_change) scale(obs_change)
## 1.975434 1.531142
## scale(aware_change):scale(nj_change) scale(nj_change):scale(obs_change)
## 1.082606 1.156170Robustness checks: Holds in univariate model
cor.test(mindfulness_changes$dr_change, mindfulness_changes$aware_change)##
## Pearson's product-moment correlation
##
## data: mindfulness_changes$dr_change and mindfulness_changes$aware_change
## t = 3.3249, df = 91, p-value = 0.001276
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1344312 0.4993550
## sample estimates:
## cor
## 0.3291264Depression change
dd1_shared <- dd1 %>% filter(ID %in% dd2$ID)
dd2_shared <- dd2 %>% filter(ID %in% dd1$ID)
if (all(dd1_shared$ID==dd2_shared$ID)) {
depr_changes <- data.frame(dd1_shared$ID,
dd2_shared$dass_sum_score-dd1_shared$dass_sum_score) %>%
setNames(c("ID", "depr_change"))
}
hist(depr_changes$depr_change)
Subset map ests…
d1me_shared_2 <- d1me %>% filter(ID %in% depr_changes$ID)… and endorsement changes
diff_ch_subs_d <- pends_p %>% filter(subj_idx %in% depr_changes$ID)Depression change
# Add drift rate change regressor
if (all(d1me_shared_2$ID==depr_changes$ID)) {
depr_changes$dr_change <- d1me_shared_2$v_val_ctr.sess_ctr_
}# Add drift rate change regressor
if (all(diff_ch_subs_d$subj_idx==depr_changes$ID)) {
depr_changes$ends_change <- diff_ch_subs_d$diff_in_diff_pn_ends
}cor.test(depr_changes$dr_change, depr_changes$depr_change)##
## Pearson's product-moment correlation
##
## data: depr_changes$dr_change and depr_changes$depr_change
## t = -2.0458, df = 91, p-value = 0.04367
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.396459540 -0.006244396
## sample estimates:
## cor
## -0.2096859cor.test(depr_changes$ends_change, depr_changes$depr_change)##
## Pearson's product-moment correlation
##
## data: depr_changes$ends_change and depr_changes$depr_change
## t = -1.8054, df = 91, p-value = 0.07431
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.37544557 0.01844734
## sample estimates:
## cor
## -0.1859598assign_short_2 <- extra_df_full %>% filter(ID %in% depr_changes$ID)testit::assert(assign_short_2$ID == depr_changes$ID)
depr_changes$tx_cond <- assign_short_2$Treatmentdepr_changes$tx_cond_f <- factor(depr_changes$tx_cond, levels=c("OM", "FA", "MBCT"))
contrasts(depr_changes$tx_cond_f) <- contr.sum(3)Check dependence on treatment condition
summary(m1 <- lm(dr_change ~ depr_change, data=depr_changes))##
## Call:
## lm(formula = dr_change ~ depr_change, data = depr_changes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.56935 -0.12910 -0.00318 0.12334 0.41858
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.090890 0.025520 3.561 0.000589 ***
## depr_change -0.006563 0.003208 -2.046 0.043667 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1867 on 91 degrees of freedom
## Multiple R-squared: 0.04397, Adjusted R-squared: 0.03346
## F-statistic: 4.185 on 1 and 91 DF, p-value: 0.04367summary(m1 <- lm(dr_change ~ tx_cond_f/depr_change, data=depr_changes))##
## Call:
## lm(formula = dr_change ~ tx_cond_f/depr_change, data = depr_changes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39715 -0.10945 -0.00604 0.10149 0.41000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.096127 0.025238 3.809 0.00026 ***
## tx_cond_f1 -0.026327 0.037552 -0.701 0.48513
## tx_cond_f2 -0.006310 0.034210 -0.184 0.85410
## tx_cond_fOM:depr_change 0.002326 0.006300 0.369 0.71287
## tx_cond_fFA:depr_change -0.011286 0.006154 -1.834 0.07007 .
## tx_cond_fMBCT:depr_change -0.008598 0.004496 -1.912 0.05912 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1808 on 87 degrees of freedom
## Multiple R-squared: 0.143, Adjusted R-squared: 0.09377
## F-statistic: 2.904 on 5 and 87 DF, p-value: 0.01796dc_vvs <- ggplot(depr_changes, aes(dr_change, depr_change)) +
geom_point(size=8, pch=21, alpha=.5, fill="gray57") +
#ggtitle(paste0( "r=-.30 p < .001 ")) +
#geom_smooth(method="lm", color="gray57") +
ga + ap +
geom_smooth(method="lm", color="gray57", alpha=.5) +
ga + theme(plot.title = element_text(margin=margin(b=-3),
size=30, hjust=.15, vjust=.6)) + theme(axis.title = element_text(size=35), axis.text=element_text(size=35)) +
ylab(" depression \n (via DASS)") + xlab("positive > negative \n drift rate * time estimates") +
theme(plot.title = element_text(margin=margin(b=5),
size=27, hjust=.15, vjust=.6)) + xlim(-.45, .75) +
geom_text(x=.55, y=0, size=15, label="*") #+ ylim(0, 37)
dc_vvs## `geom_smooth()` using formula 'y ~ x'
#ggsave("../paper/figs/depr-change_v-val_sess.png", dc_vvs, width=9, height=7)Robustness check: Adding depression as covariate in mindfulness model
if (all(mindfulness_changes$ID==depr_changes$ID)) {
mindfulness_changes <- data.frame(mindfulness_changes, depr_changes$depr_change)
}
summary(
change_2m <-
lm(dr_change ~
scale(aware_change)*scale(nj_change) +
scale(nr_change) + scale(depr_changes.depr_change),
data=mindfulness_changes)
)##
## Call:
## lm(formula = dr_change ~ scale(aware_change) * scale(nj_change) +
## scale(nr_change) + scale(depr_changes.depr_change), data = mindfulness_changes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53744 -0.11422 -0.00769 0.11486 0.44203
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.128345 0.020440 6.279 1.3e-08 ***
## scale(aware_change) 0.057527 0.022503 2.556 0.0123 *
## scale(nj_change) 0.019432 0.024566 0.791 0.4311
## scale(nr_change) -0.028330 0.022014 -1.287 0.2015
## scale(depr_changes.depr_change) -0.016340 0.021043 -0.777 0.4395
## scale(aware_change):scale(nj_change) -0.006945 0.016235 -0.428 0.6699
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1812 on 87 degrees of freedom
## Multiple R-squared: 0.139, Adjusted R-squared: 0.08951
## F-statistic: 2.809 on 5 and 87 DF, p-value: 0.02123cor.test(mindfulness_changes$depr_changes.depr_change, mindfulness_changes$aware_change)##
## Pearson's product-moment correlation
##
## data: mindfulness_changes$depr_changes.depr_change and mindfulness_changes$aware_change
## t = -3.1211, df = 91, p-value = 0.002415
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4840102 -0.1145071
## sample estimates:
## cor
## -0.3109618Relationship between change in different mindfulness facets and depression change
cor.test(mindfulness_changes$aware_change, mindfulness_changes$depr_changes.depr_change)##
## Pearson's product-moment correlation
##
## data: mindfulness_changes$aware_change and mindfulness_changes$depr_changes.depr_change
## t = -3.1211, df = 91, p-value = 0.002415
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4840102 -0.1145071
## sample estimates:
## cor
## -0.3109618cor.test(mindfulness_changes$nj_change, mindfulness_changes$depr_changes.depr_change)##
## Pearson's product-moment correlation
##
## data: mindfulness_changes$nj_change and mindfulness_changes$depr_changes.depr_change
## t = -4.1273, df = 91, p-value = 8.117e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5558350 -0.2103992
## sample estimates:
## cor
## -0.3970883cor.test(mindfulness_changes$nr_change, mindfulness_changes$depr_changes.depr_change)##
## Pearson's product-moment correlation
##
## data: mindfulness_changes$nr_change and mindfulness_changes$depr_changes.depr_change
## t = -0.82196, df = 91, p-value = 0.4132
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2845783 0.1199595
## sample estimates:
## cor
## -0.08584661Robustness check: Depression and DDM – still relevant?
# cor.test(mindfulness_changes$dr_change_ddm,
# mindfulness_changes$depr_changes.depr_change)Results: Split-half and test-retest reliability of model parameters
These are also too big to put in Public_Data
even_ddm <- read.csv("../../model_res/traces/DDM_split_half_even_wtrialwise_poutlier057074.csv")
odd_ddm <- read.csv("../../model_res/traces/DDM_split_half_odd__wtrialwise_poutlier052866.csv")
pre_t_d <- read.csv("./../../model_res/traces/DDM_test_retest_PRE_s-val_wtrialwise_poutlier059850.csv")
post_t_d <- read.csv("./../../model_res/traces/DDM_test_retest_POST_s-val_wtrialwise_poutlier055680.csv")even_emp <- read.csv("../../data/cleaned_files/s_bdf_even.csv")
odd_emp <- read.csv("../../data/cleaned_files/s_bdf_odd.csv")
# even_ddm <- read.csv("../../model_res/traces/DDM_split_half_even2558.csv")
# odd_ddm <- read.csv("../../model_res/traces/DDM_split_half_odd3780.csv")GetICCData <- function(vec1, vec2) {
out <- psych::ICC(data.table(vec1, vec2))
data.table(out$results$ICC[2], out$results$`lower bound`[2], out$results$`upper bound`[2]) %>%
setNames(c("ICC", "95_ci_lb", "95_ci_ub"))
}Behavioral effects
Split-half
even_rt <- even_emp %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)
odd_rt <- odd_emp %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)
even_rt_n <- even_emp %>% filter(valence=="negative") %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)
odd_rt_n <- odd_emp %>% filter(valence=="negative") %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)
even_rt_p <- even_emp %>% filter(valence=="positive") %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)
odd_rt_p <- odd_emp %>% filter(valence=="positive") %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)Calculate valence rt diffs
if (all(even_rt_n$subj_idx==odd_rt_n$subj_idx)) {
even_v_rt_diff <- even_rt_p$rt - even_rt_n$rt
odd_v_rt_diff <- odd_rt_p$rt - odd_rt_n$rt
cor.test(even_v_rt_diff, odd_v_rt_diff)
psych::ICC(data.table(even_v_rt_diff, odd_v_rt_diff))
}## boundary (singular) fit: see help('isSingular')## Call: psych::ICC(x = data.table(even_v_rt_diff, odd_v_rt_diff))
##
## Intraclass correlation coefficients
## type ICC F df1 df2 p lower bound upper bound
## Single_raters_absolute ICC1 0.36 2.1 96 97 0.00013 0.21 0.50
## Single_random_raters ICC2 0.36 2.1 96 96 0.00013 0.21 0.50
## Single_fixed_raters ICC3 0.36 2.1 96 96 0.00013 0.21 0.50
## Average_raters_absolute ICC1k 0.53 2.1 96 97 0.00013 0.34 0.66
## Average_random_raters ICC2k 0.53 2.1 96 96 0.00013 0.34 0.66
## Average_fixed_raters ICC3k 0.53 2.1 96 96 0.00013 0.34 0.66
##
## Number of subjects = 97 Number of Judges = 2
## See the help file for a discussion of the other 4 McGraw and Wong estimates,even_pends <- even_emp %>% filter( valence=="positive") %>%
group_by(subj_idx) %>% summarize_at(c("response"), mean, na.rm=TRUE)
even_nends <- even_emp %>% filter(valence=="negative") %>%
group_by(subj_idx) %>% summarize_at(c("response"), mean, na.rm=TRUE)
even_pos_min_neg <- even_pends$response - even_nends$response
odd_pends <- odd_emp %>% filter(valence=="positive") %>%
group_by(subj_idx) %>% summarize_at(c("response"), mean, na.rm=TRUE)
odd_nends <- odd_emp %>% filter(valence=="negative") %>%
group_by(subj_idx) %>% summarize_at(c("response"), mean, na.rm=TRUE)
odd_pos_min_neg <- odd_pends$response - odd_nends$responseTest-retest
pre_rt <- bx_df %>% filter(session=="Pre") %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)
post_rt <- bx_df %>% filter(session=="Post") %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)
pre_rt_n <- bx_df %>% filter(session=="Pre", valence=="negative") %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)
post_rt_n <- bx_df %>% filter(session=="Post", valence=="negative") %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)
pre_rt_p <- bx_df %>% filter(session=="Pre", valence=="positive") %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)
post_rt_p <- bx_df %>% filter(session=="Post", valence=="positive") %>%
group_by(subj_idx) %>% summarize_at(c("rt"), mean, na.rm=TRUE)if (all(pre_rt$subj_idx == post_rt$subj_idx)) cor.test(pre_rt$rt, post_rt$rt); psych::ICC(data.table(pre_rt$rt, post_rt$rt))##
## Pearson's product-moment correlation
##
## data: pre_rt$rt and post_rt$rt
## t = 8.6151, df = 94, p-value = 1.611e-13
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5349965 0.7631006
## sample estimates:
## cor
## 0.6642369## Call: psych::ICC(x = data.table(pre_rt$rt, post_rt$rt))
##
## Intraclass correlation coefficients
## type ICC F df1 df2 p lower bound upper bound
## Single_raters_absolute ICC1 0.63 4.4 95 96 2.4e-12 0.51 0.72
## Single_random_raters ICC2 0.64 4.9 95 95 7.0e-14 0.50 0.73
## Single_fixed_raters ICC3 0.66 4.9 95 95 7.0e-14 0.56 0.75
## Average_raters_absolute ICC1k 0.77 4.4 95 96 2.4e-12 0.68 0.84
## Average_random_raters ICC2k 0.78 4.9 95 95 7.0e-14 0.67 0.85
## Average_fixed_raters ICC3k 0.80 4.9 95 95 7.0e-14 0.72 0.86
##
## Number of subjects = 96 Number of Judges = 2
## See the help file for a discussion of the other 4 McGraw and Wong estimates,if (all(pre_rt_n$subj_idx == post_rt_n$subj_idx)) cor.test(pre_rt_n$rt, post_rt_n$rt); psych::ICC(data.table(pre_rt_n$rt, post_rt_n$rt)); test <- "ok"##
## Pearson's product-moment correlation
##
## data: pre_rt_n$rt and post_rt_n$rt
## t = 7.7012, df = 94, p-value = 1.346e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4815411 0.7312768
## sample estimates:
## cor
## 0.6219797## Call: psych::ICC(x = data.table(pre_rt_n$rt, post_rt_n$rt))
##
## Intraclass correlation coefficients
## type ICC F df1 df2 p lower bound upper bound
## Single_raters_absolute ICC1 0.57 3.7 95 96 3.2e-10 0.45 0.68
## Single_random_raters ICC2 0.59 4.2 95 95 7.0e-12 0.44 0.70
## Single_fixed_raters ICC3 0.62 4.2 95 95 7.0e-12 0.50 0.71
## Average_raters_absolute ICC1k 0.73 3.7 95 96 3.2e-10 0.62 0.81
## Average_random_raters ICC2k 0.74 4.2 95 95 7.0e-12 0.61 0.82
## Average_fixed_raters ICC3k 0.76 4.2 95 95 7.0e-12 0.67 0.83
##
## Number of subjects = 96 Number of Judges = 2
## See the help file for a discussion of the other 4 McGraw and Wong estimates,if (all(pre_rt_p$subj_idx == post_rt_p$subj_idx)) cor.test(pre_rt_p$rt, post_rt_p$rt); psych::ICC(data.table(pre_rt_p$rt, post_rt_p$rt)); test2 <- "ok"##
## Pearson's product-moment correlation
##
## data: pre_rt_p$rt and post_rt_p$rt
## t = 8.628, df = 94, p-value = 1.513e-13
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5357038 0.7635144
## sample estimates:
## cor
## 0.6647905## Call: psych::ICC(x = data.table(pre_rt_p$rt, post_rt_p$rt))
##
## Intraclass correlation coefficients
## type ICC F df1 df2 p lower bound upper bound
## Single_raters_absolute ICC1 0.64 4.6 95 96 5.4e-13 0.53 0.73
## Single_random_raters ICC2 0.65 4.9 95 95 6.9e-14 0.53 0.74
## Single_fixed_raters ICC3 0.66 4.9 95 95 6.9e-14 0.56 0.75
## Average_raters_absolute ICC1k 0.78 4.6 95 96 5.4e-13 0.69 0.84
## Average_random_raters ICC2k 0.79 4.9 95 95 6.9e-14 0.69 0.85
## Average_fixed_raters ICC3k 0.80 4.9 95 95 6.9e-14 0.72 0.86
##
## Number of subjects = 96 Number of Judges = 2
## See the help file for a discussion of the other 4 McGraw and Wong estimates,Calculate valence rt diffs. A bit worse but even it’s not bad
if (test=="ok" && test2=="ok") {
pre_v_rt_diff <- pre_rt_p$rt - pre_rt_n$rt
post_v_rt_diff <- post_rt_p$rt - post_rt_n$rt
cor.test(pre_v_rt_diff, post_v_rt_diff)
psych::ICC(data.table(pre_v_rt_diff, post_v_rt_diff))
} else {
cat("Error with data alignment—see above cell.")
}## Call: psych::ICC(x = data.table(pre_v_rt_diff, post_v_rt_diff))
##
## Intraclass correlation coefficients
## type ICC F df1 df2 p lower bound upper bound
## Single_raters_absolute ICC1 0.42 2.5 95 96 7.7e-06 0.27 0.55
## Single_random_raters ICC2 0.42 2.5 95 95 7.2e-06 0.28 0.55
## Single_fixed_raters ICC3 0.42 2.5 95 95 7.2e-06 0.28 0.55
## Average_raters_absolute ICC1k 0.59 2.5 95 96 7.7e-06 0.43 0.71
## Average_random_raters ICC2k 0.59 2.5 95 95 7.2e-06 0.43 0.71
## Average_fixed_raters ICC3k 0.60 2.5 95 95 7.2e-06 0.43 0.71
##
## Number of subjects = 96 Number of Judges = 2
## See the help file for a discussion of the other 4 McGraw and Wong estimates,pre_pends <- bx_df %>% filter(session=="Pre", valence=="positive") %>%
group_by(subj_idx) %>% summarize_at(c("response"), mean, na.rm=TRUE)
pre_nends <- bx_df %>% filter(session=="Pre", valence=="negative") %>%
group_by(subj_idx) %>% summarize_at(c("response"), mean, na.rm=TRUE)
pre_pos_min_neg <- pre_pends$response - pre_nends$response
post_pends <- bx_df %>% filter(session=="Post", valence=="positive") %>%
group_by(subj_idx) %>% summarize_at(c("response"), mean, na.rm=TRUE)
post_nends <- bx_df %>% filter(session=="Post", valence=="negative") %>%
group_by(subj_idx) %>% summarize_at(c("response"), mean, na.rm=TRUE)
post_pos_min_neg <- post_pends$response - post_nends$response
# post_pends <- bx_df %>% filter(session=="Post", valence=="positive") %>%
# group_by(subj_idx) %>% summarize_at(c("response"), mean, na.rm=TRUE)
#
# pre_pends <- bx_df %>% filter(session=="Pre", valence=="positive") %>%
# group_by(subj_idx) %>% summarize_at(c("response"), mean, na.rm=TRUE)
#
# post_pends <- bx_df %>% filter(session=="Post", valence=="positive") %>%
# group_by(subj_idx) %>% summarize_at(c("response"), mean, na.rm=TRUE)Following Hedge et al. 17 using ICC2
Reliabilities were calculated using Intraclass Correlation Coefficients (ICC) using a two-way random effects model for absolute agreement. In the commonly cited Shrout and Fleiss (1979; see also McGraw & Wong, 1996) nomenclature, this corresponds to ICC (2,1). This form of the ICC is sensitive to differences between session means. In Supplementary Material A, we perform further analyses to account for potential outliers and distributional assumptions. The choice of statistic does not affect our conclusions. We report reliabilities separately for Studies 1 and 2 in the main text so that consistency across samples can be observed. We combine the studies in supplementary analyses.
Split half
mf_ICCs_sh <- rbind(
data.frame(GetICCData(even_rt$rt, odd_rt$rt), "mean \nreaction time") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(even_v_rt_diff, odd_v_rt_diff), "difference in \nreaction time by valence") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(even_nends$response, odd_nends$response), "proportion neg. endorsements") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(even_pends$response, odd_pends$response), "proportion pos. endorsements") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(even_pos_min_neg, odd_pos_min_neg), "ratio of positive:negative \n endorsements") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var"))
)## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')# mf_ICCs_sh$lb <- mf_ICCs_sh$ICC-mf_ICCs_sh$ci_lb
# mf_ICCs_sh$ub <- mf_ICCs_sh$ci_ub - mf_ICCs_sh$ICC
mf_ICCs_sh$type <- "split_half"mf_ICCs_sh## ICC ci_lb ci_ub var
## 1 0.9155235 0.8837085 0.9389370 mean \nreaction time
## 2 0.3603387 0.2059824 0.4972983 difference in \nreaction time by valence
## 3 0.5495007 0.4214780 0.6561456 proportion neg. endorsements
## 4 0.4546658 0.3115902 0.5776796 proportion pos. endorsements
## 5 0.6119708 0.4958958 0.7066414 ratio of positive:negative \n endorsements
## type
## 1 split_half
## 2 split_half
## 3 split_half
## 4 split_half
## 5 split_halfTest retest
mf_ICCs_tr <- rbind(
data.frame(GetICCData(pre_rt$rt, post_rt$rt), "mean \nreaction time") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(pre_v_rt_diff, post_v_rt_diff), "difference in \nreaction time by valence") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(pre_nends$response, post_nends$response), "proportion neg. endorsements") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(pre_pends$response, post_pends$response), "proportion pos. endorsements") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(pre_pos_min_neg, post_pos_min_neg), "ratio of positive:negative \n endorsements") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var"))
)
# mf_ICCs_tr$lb <- mf_ICCs_tr$ICC-mf_ICCs_tr$ci_lb
# mf_ICCs_tr$ub <- mf_ICCs_tr$ci_ub - mf_ICCs_tr$ICC
mf_ICCs_tr$type <- "test_retest"mf_ICCs_tr## ICC ci_lb ci_ub var
## 1 0.6356479 0.5028770 0.7346966 mean \nreaction time
## 2 0.4233395 0.2758852 0.5515546 difference in \nreaction time by valence
## 3 0.5412898 0.3742585 0.6660359 proportion neg. endorsements
## 4 0.5640652 0.4313549 0.6714392 proportion pos. endorsements
## 5 0.5352843 0.3614335 0.6635753 ratio of positive:negative \n endorsements
## type
## 1 test_retest
## 2 test_retest
## 3 test_retest
## 4 test_retest
## 5 test_retestmf_iccs <- rbind(
data.frame(mf_ICCs_sh, "label"="split-half"),
data.frame(mf_ICCs_tr, "label"="test-retest")
)Model-based measures
Split half
DDM
d_even_me <- GetMapEsts(GetSubjTraces(even_ddm), format="short")
d_odd_me <- GetMapEsts(GetSubjTraces(odd_ddm), format="short")
# Combine the datasets
if (all(d_even_me$ID == d_odd_me$ID)) {
c_me_d <- data.table(d_even_me, d_odd_me[2:ncol(d_odd_me)] %>%
setNames(paste0(names(d_odd_me[2:ncol(d_odd_me)]), "odd")))
}
mb_d_ICCs_sh <- rbind(
data.frame(GetICCData(c_me_d$a_, c_me_d$a_odd), "threshold") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(c_me_d$t_, c_me_d$t_odd), "non-decision time") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
# Dummy var for this since not in DDM
# The one that generates convergence error.. also has lowest ICC
data.frame(GetICCData(InvLogit(c_me_d$z_), InvLogit(c_me_d$z_odd)), "starting point bias") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(c_me_d$v_Intercept_, c_me_d$v_Intercept_odd), "drift rate - intercept") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(c_me_d$v_val_ctr_, c_me_d$v_val_ctr_odd), "drift rate - pos. vs. neg. \nregressor") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var"))
#data.frame(GetICCData(data.table(even_rt$rt, odd_rt$rt), "RT mean") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var"))
)## boundary (singular) fit: see help('isSingular')# mb_d_ICCs_sh$lb <- mb_d_ICCs_sh$ICC-mb_d_ICCs_sh$ci_lb
# mb_d_ICCs_sh$ub <- mb_d_ICCs_sh$ci_ub-mb_d_ICCs_sh$ICCmb_d_ICCs_sh## ICC ci_lb ci_ub var
## 1 0.647664393 0.536589071 0.73616986 threshold
## 2 0.777602330 0.701923310 0.83596992 non-decision time
## 3 0.003412813 -0.002956936 0.01331431 starting point bias
## 4 0.396511669 0.220357617 0.54120637 drift rate - intercept
## 5 0.671413324 0.564663987 0.75534506 drift rate - pos. vs. neg. \nregressorTest retest
DDM
#post_t_d <- read.csv("./../model_res/traces/DDM_test_retest_POST_m0_sret_v_val2062.csv") # Sub in the correct one
d_pre_me <- GetMapEsts(GetSubjTraces(pre_t_d), format="short")
d_post_me <- GetMapEsts(GetSubjTraces(post_t_d), format="short")
# Combine the datasets
if (all(d_pre_me$ID == d_post_me$ID)) {
c_me_d <- data.table(d_pre_me, d_post_me[2:ncol(d_post_me)] %>%
setNames(paste0(names(d_post_me[2:ncol(d_post_me)]), "post")))
}
mb_d_ICCs <- rbind(
data.frame(GetICCData(c_me_d$a_, c_me_d$a_post), "threshold") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(c_me_d$t_, c_me_d$t_post), "non-decision time") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
# Dummy var for this since not in DDM
# The one that generates convergence error.. also has lowest ICC
data.frame(GetICCData(InvLogit(c_me_d$z_), InvLogit(c_me_d$z_post)), "starting point bias") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(c_me_d$v_Intercept_, c_me_d$v_Intercept_post), "drift rate - intercept") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var")),
data.frame(GetICCData(c_me_d$v_val_ctr_, c_me_d$v_val_ctr_post), "drift rate - pos. vs. neg. \nregressor") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var"))
#data.frame(GetICCData(data.table(pre_rt$rt, post_rt$rt), "RT mean") %>% setNames(c("ICC", "ci_lb", "ci_ub", "var"))
)## boundary (singular) fit: see help('isSingular')# mb_d_ICCs$lb <- mb_d_ICCs$ICC-mb_d_ICCs$ci_lb
# mb_d_ICCs$ub <- mb_d_ICCs$ci_ub-mb_d_ICCs$ICCddm_iccs <- rbind(
data.frame(mb_d_ICCs_sh, "label"="split-half"),
data.frame(mb_d_ICCs, "label"="test-retest")
)mf_p <- ggplot(mf_iccs, aes(x=var, y=ICC, group=label, fill=label)) +
geom_errorbar(aes(ymin=ci_lb, ymax=ci_ub), position=position_dodge(width=.4), size=2, width=0) +
geom_hline(yintercept=c(0, .25, .5, .75, 1)) +
geom_point(alpha=.9, alpha=.9, size=8, pch=21, position=position_dodge(width=.4)) +
#geom_hline(yintercept = 0, linetype="dotted") +
theme(axis.text.x = element_text(angle=25, hjust=1, size=16)) +
ga + ap + ylab("") + ylim(-.1, 1) + xlab("") +
#ggtitle("Model-derived measures") +
#theme(axis.text.y = element_blank(), axis.ticks.y=element_blank()) +
theme(plot.title = element_text(margin=margin(b=-25),
size=30, hjust=.15, vjust=.6)) +
scale_fill_manual(values=c("grey24", "gray65")) + lp #+## Warning: Duplicated aesthetics after name standardisation: alpha # ggtitle("Behavioral measures") +
# theme(plot.title = element_text(margin=margin(b=-25),
# size=30, hjust=.15, vjust=.6))
#mf_pmb_p <- ggplot(ddm_iccs, aes(x=var, y=ICC, group=label, fill=label)) +
geom_errorbar(aes(ymin=ci_lb, ymax=ci_ub), position=position_dodge(width=.4), size=2, width=0) +
geom_hline(yintercept=c(0, .25, .5, .75, 1)) +
geom_point(alpha=.9, alpha=.9, size=8, pch=21, position=position_dodge(width=.4)) +
#geom_hline(yintercept = 0, linetype="dotted") +
theme(axis.text.x = element_text(angle=25, hjust=1, size=16)) +
ga + ap + ylab("") + ylim(-.1, 1) + xlab("") +
#ggtitle("Model-derived measures") +
#theme(axis.text.y = element_blank(), axis.ticks.y=element_blank()) +
theme(plot.title = element_text(margin=margin(b=-25),
size=30, hjust=.15, vjust=.6)) +
scale_fill_manual(values=c("tan1", "chocolate")) + lp## Warning: Duplicated aesthetics after name standardisation: alphaiccs_p <- mf_p / mb_p +
plot_annotation(
title = 'Stability of behavioral (top) and model-based (bottom) measures',
subtitle='ICC point estimates and 95% CIs',
theme = theme(plot.title = element_text(size = 25, hjust=.5),
plot.subtitle = element_text(size = 20, hjust=.05))
)
iccs_p
#ggsave("../paper/figs/icc_plot.png", iccs_p, width=11, height=10)3.25 - Note that after R upgrade factor rule changes must have led to change in order on x-axis compared to version in figure in paper, but spot-checked to confirm all results were the same.